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Preface 


The  purpose  of  this  study  was  to  explore  the  feasibility 
of  using  the  finite  element  method  in  a  numerical  solution  of 
the  transport  equation.  For  problems  in  science  and  engi¬ 
neering,  the  finite  element  method  can  offer  significant 
advantages  over  competing  numerical  techniques,  such  as  the 
method  of  finite  difference.  In  general,  it  requires  that 
the  differential  equation  modeling  the  physical  system  be 
mathematically  self-adjoint.  The  finite  element  method  has 
not  been  used  in  the  past  on  the  transport  equation  for  pre¬ 
cisely  this  reason:  it  is  not,  in  its  current  formulation,  a 
self-adjoint  equation.  In  this  study,  the  transport  equation 
was  recast  into  a  form  which,  while  more  complex,  is  indeed 
self-adjoint.  The  finite  element  equations  were  derived  for 
the  one-dimensional  case  with  isotropic  scatter.  This  proved 
to  be  one  of  the  challenges  in  this  study  due  to  the  non¬ 
local  character  of  the  transport  equation.  The  term  "non- 
local"  means  that  the  solution  at  a  given  point  is  dependent 
upon  points  which,  in  phase  space,  are  far  removed  from  the 
point  of  interest.  Lastly,  a  numerical  solution  was  written 
in  FORTRAN  and  implemented  on  a  VAX  11/780.  The  results  were 
compared  against  a  benchmark  solution  using  a  recently  devel¬ 
oped  solution  technique  called  Ln. 


As  in  all  research  efforts,  several  people  played  impor¬ 
tant  roles.  I  would  like  to  particularly  thank  Dr.  Donn 
Shankland  for  the  technical  guidance  he  provided.  His  con¬ 
tribution  to  this  research  effort  was  significant.  I  would 
also  like  to  express  my  gratitude  to  Dr.  John  Jones  for  his 
suggestions  on  the  finite  element  method.  Finally,  I  would 
like  to  thank  my  wife,  Jamie,  whose  love  and  understanding 
endured  through  this  graduate  program. 
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Abstract 


A  self-adjoint  form  of  the  transport  equation  was 
derived,  and  expressed  as  an  extremization  integral.  The 
finite  element  equations  were  derived  from  the  extremization 
integral  for  the  one-dimensional  time  independent  homogeneous 
transport  equation  with  isotropic  scatter.  These  equations 
were  implemented  in  FORTRAN  on  a  VAX  11/780  and  used  to  solve 
a  simple  benchmark  problem. 

The  finite  element  solution,  for  a  small  mesh  of  32  ele¬ 
ments,  was  compared  to  results  from  a  numerical  technique 
known  as  Ln.  The  solutions  differed  by  about  35  percent. 
Larger  meshes  were  not  run  because  an  automatic  mesh  refine¬ 
ment  routine  was  not  available.  The  large  difference  between 
the  Ln  solution  and  the  finite  difference  solution  is  attri¬ 
buted  to  residual  errors  in  the  coding  of  the  finite  element 


equations  for  the  scatter  terms. 


FINITE  ELEMENT  SOLUTION  OF  A  SELF-ADJOINT 
TRANSPORT  EQUATION  IN  ONE  DIMENSION 


I .  Introduction 

« 

Background 

The  transport  equation  is  a  remarkably  successful  des¬ 
cription  of  transport  phenomena.  It  has  been  used  on  a  wide 
range  of  problems  in  such  diverse  fields  as  nuclear  reactors, 
astrophysics,  charged  particle  transport,  rarefied  gas  dynam¬ 
ics,  and  even  traffic  flow. 

A  variety  of  solution  techniques  have  been  applied  to 
the  transport  equation  (4:2),  which  may  be  roughly  divided 
into  four  classes:  analytic,  approximation,  numerical,  and 
simulation.  Except  for  the  very  simplest  problems,  the 
analytic  solutions  tend  to  be  specialized  to  some  particular 
field,  or  problem,  and  therefore  lack  generality.  The  stan¬ 
dard  approximation  schemes  are  complicated  by  the  non-self- 
adjoint  nature  of  the  transport  equation,  and  they  often 
require  considerable  modification.  Simulations,  such  as 
Monte  Carlo,  are  capable  of  very  high  accuracy,  but  the  cost 
is  often  prohibitive  due  to  the  large  number  of  computer 
operations  required. 

Numerical  techniques  offer  a  good  engineering  compro¬ 
mise.  They  are  reasonably  accurate  at  a  reasonable  cost  and 


are  generally  applicable  over  a  wide  range  of  problems. 

The  finite  element  method  is  a  relatively  recent  (7:vii) 
numerical  technique  which  has  been  successfully  applied  in 
many  fields.  It  has  been  applied  to  the  transport  equation, 
without  notable  success,  due  to  the  non-self-adjoint  nature 
of  the  transport  equation  (4:479-504). 


Objective 

The  objective  of  this  study  was  to  recast  the  transport 
equation  into  a  form  which  is  self-adjoint,  and  then  apply 
the  finite  element  method  toward  its  solution. 


Scope 

The  transport  equation  was  recast  into  a  self-adjoint 
form  and  expressed  as  an  extremization  principle.  The  finite 
element  equations  were  derived  from  this  extremization  prin¬ 
ciple  for  the  homogeneous  transport  equation  with  isotropic 
scatter.  They  were  then  implemented  in  FORTRAN  on  a  VAX 
11/780  for  a  simple  benchmark  problem.  Two  sets  of  computer 
runs  were  done.  In  the  first  set,  the  results  for  the  trans¬ 
port  equation  without  scatter  were  compared  to  an  analytic 
solution.  In  the  second  set,  the  results  for  the  transport 
equation  with  scatter  were  compared  to  the  results  from 
another  numerical  technique  called  Ln. 


Assumptions 


The  assumptions  in  this  study  are  those  commonly  used 


in  deriving  the  transport  equation,  those  used  in  making  the 
finite  element  approximation,  and  those  used  in  defining  the 
benchmark  problem.  They  are  described  in  the  relevant  sec¬ 
tions  of  this  thesis. 

Sequence  of  Presentation 

The  next  chapter  reviews  the  transport  equation  and  pre¬ 
sents  the  derivation  of  the  self-adjoint  form.  In  Chapter 
III,  the  finite  element  equations  are  derived.  The  computer 
implementation  is  described  in  Chapter  IV.  Chapter  V  pre¬ 
sents  the  results  of  the  computer  program  and  Chapter  VI 
gives  the  conclusions  and  recommendations. 
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The  term  "transport  theory"  is  commonly  used  to  refer  to 
the  mathematical  description  of  the  motion  of  particles 
through  a  host  medium.  It  differs  from  the  usual  approaches 
used  in  classical  physics  because  it  is  a  particle,  rather 
than  a  continuum,  theory  of  matter.  The  concept  of  a  con¬ 
tinuous  field  is,  however,  employed  in  transport  theory  in 
probability  fields  for  some  of  the  properties  of  the  system 
of  particles. 

There  are  several  common  assumptions  in  transport 
theory.  In  this  development,  the  following  assumptions  were 
made: 

The  particles  interact  only  with  the  host  medium. 

There  is  no  correlation  between  particles. 

Only  the  particle's  position  and  velocity  are  sig¬ 
nificant. 

All  interactions  with  the  host  medium  occur  at  a 
point  well  localized  in  space  with  respect  to  the 
mean  free  path  and  well  localized  in  time  with 
respect  to  the  mean  time  of  flight. 

Transport  theory  rests  on  several  key  concepts  which  are 
reviewed  briefly  here.  Greater  detail  can  be  found  in 
Transport  Theory,  by  Duderstadt  and  Martin  (4:1-19). 

The  first  and  most  central  concept  is  that  of  a  proba¬ 
bility  field  which  describes  the  expected  density  of 


particles  in  configuration  space.  Specifically,  the  particle 
density  is  defined  by 

ML  ?,t)Sr  —  expected  number  of  particles  in 
about  at  time  t. 


Similarly,  the  phase  space  particle  density  is  defined  by 


time  t. 


The  two  are  related  by 

M(r]t)  =  / y,  t)  J?\r  ( 1 ) 

The  second  key  concept  is  the  representation  of  inter¬ 
actions  between  the  transporting  particles  and  the  host 
medium.  Two  terms  are  used  in  this  concept.  They  are  "mean 
free  path"  (mfp)  and  "macroscopic  cross-section."  They  are 
inverses  of  each  other  and  are  defined  by 


probability  of  particle  interaction 
per  unit  distance  traveled  by  a 
particle  of  velocity  Xr  at  position 

/  . 


As  a  further  refinement,  the  collision  kernel  is  defined  as 
the 


probability  per  unit  distance  trav¬ 
eled  that^  an  incident  particle  of 
velocity  V  will  collide  producing 
a  secondary  particle  at  velocity  \r  . 


(Vp  )'=£(?£;* 


Mathematically 


— 7/  — >  \ 

V-j 


(2) 


where  c  is  the  mean  number  of  secondary  particles  emitted  per 
collision  event  and  is  defined  by 


tlr.v  J  S  mean  number  of  secondary  particles 
emitted  in  a  collision  event  experi¬ 
enced  by  an  incident  particle  with 
velocity  *[?  at  position  lr. 


and  where  f  is  the  scattering  probability  function  defined  by 

probability  that  any  secondary  par¬ 
ticles  induced  by  an  incident  parti¬ 
cle  with  velocity  v*  will  be  emitted 
with  velocity  Xr  in  <dL*\r. 

Note  that 

ff  -  I  (3) 

When  the  macroscopic  cross-section  is  independent  of  position 
and  velocity  it  is  symbolized  by  for  total. 

In  addition  to  these  key  concepts,  there  are  a  few  sup¬ 
porting  concepts.  The  particle  scalar  flux  is  defined  by 

<P(tl  tr)  =  y-tt)£v-  (4) 

where  (ft  Cr,\s,t)  is  the  particle  scalar  phase  space  flux  and 


{ «A  a 


is  defined  by 


The  term  "flux"  is  usually  used  by  itself,  with  the  adjec¬ 
tives  dropped,  relying  on  context  to  identify  the  desired 
meaning.  This  use  of  the  term  is,  except  in  the  nuclear  com¬ 
munity,  non-standard  and  should  not  be  confused  with  the  vec¬ 
tor  flux.  The  phase  space  current  density  is  defined  as 

J( K  *)  (6) 


where 


jlTtir,i)'  di£  <£ir  s 


expected  number  of  particles  that 
cross  an  area  cLS  per  second  with 
velocity  Tr  in  et^lr  at  time  t. 


The  last  supporting  concept  is  of  a  unit  vector  in  the 
direction  of  the  velocity  vector,  called  the  cosine  vector. 
It  is  defined  by 


m  yj£j  —  t*SL<P  f  fly  j-  (7) 

where  &  and  have  their  usual  meanings  in  spherical  coordi- 

-i  -f 

nates  and  and  are  the  unit  vectors  in  cartesian 

coordinates. 

The  transport  equation  can  be  derived  by  invoking  a 
balance  condition  around  the  conservation  of  particles.  In 
an  arbitrarily  small,  fixed  region  in  space  with  volume  V, 
surface  area  S  and  number  of  particles  n,  the  balance  condi¬ 
tion  is 


ftime  rate^  (change  due^  C  change  ^  f 

<of  change?  =  <to  leakage y  +  V  due  to  r  +4s 

L  of  n  )  L  thru  S  /  (_collisions)  L 


(sources 


Writing  the  balance  condition  mathematically  gives 

^f[v [nit, K  A  = 

V  +JCL^A]**  +fr[sli:St)4Wr  (8> 

The  leakage  term  can  be  rewritten  using  Gauss's  law  and  using 
the  independence  between  the  space  and  velocity  variables  as 

[<J  (x,  y  (9> 

/,  [?■  j(?  5  5  (10) 


Thus,  the  balance  condition  becomes 


f{[%  *  -  (-£%,  -*)  M  ^  = 


(id 


Since  the  volume  is  arbitrary 


=  s 


(12) 


The  collision  term  can  be  represented  more  precisely  by 

(13) 

which,  under  the  assumption  of  uniform  host  medium,  becomes 

$t)un  ~  Jv  '£  (if  U  v)r\(r]  v;  t)  -lr$t  (y)  +)  (14) 

Substitution  gives 

f  2r*  VK  f-  2r^(zTJ)i  =  J>£(ir-'?v)i/,n(tfi%t)4ir' +  5  (is) 


Rewriting  it  in  terms  of  flux  and  cosine  vectors,  gives 

~fd£ ‘J*  £(f-yE/  £-**)#(£  t)jA+  S  (16) 

While  this  form  of  the  transport  equation  is  very  explicit, 
it  is  unwieldy  for  many  mathematical  manipulations.  As  a 
notational  aid,  an  operator  is  defined  such  that  the  trans¬ 
port  equation  becomes 

/<f> -S  =  o  (17) 

where,  in  this  instance,  the  operator  is 

/=  V  it  ~  (£-?£,*■»■*)  cf*  (18) 

As  in  all  differential  equations,  this  statement  of  the 
problem  is  incomplete  until  boundary  conditions  are  speci¬ 
fied. 

Depending  on  the  application,  several  additional  assump¬ 
tions  are  commonly  applied  to  simplify  the  mathematics  of  the 
transport  equation.  In  this  study,  four  additional  assump¬ 
tions  were  made. 

The  first  assumption  was  that  all  sources  are  mono- 
energetic  and  that  collision  events  result  in  either  absorp¬ 
tion  or  in  elastic  scatter.  Thus,  the  transport  equation 
becomes 

a:  5  (19) 

and  is  called  the  one  speed  model. 
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The  second  assumption  was  to  ignore  the  time  dependency, 
that  is,  only  the  steady  state-solution  is  of  interest. 
Thus,  the  transport  equation  becomes 


i-it*  -  ><#  =  5 


(20) 


The  third  assumption  was  that  the  region  of  interest  had 
planar  geometry,  i.e.,  that  it  could  be  modeled  in  one  dimen¬ 
sion.  Aligning  the  coordinate  system  with  the  z-axis  perpen¬ 
dicular  to  the  plane  of  symmetry  simplifies  the  divergence  to 

V<P  =  4i  e\ 


(21) 


Therefore 


(22) 


letting 


AA  *=  tOZ  0 


(23) 


and  changing  the  label  on  the  z-axis  to  x  yields 


AX' 


—  ^ 


(24) 


Additionally,  expanding  the  collision  kernel  yields 

)  —  t  5-rf( A->^)  =  £5  {(A  '-*J\  J  (25) 


so  the  one-dimensional  transport  equation  becomes 

M  +■  i.t<P  ~  f /(A->A)  <fU,  -  c. 


(26) 


r';/. 


The  fourth  and  final  assumption  was  that  of  isotropic 


scatter.  Recall  that 


ffisi’-ysi)  eC*  —  1 


(27) 


since 


—  5  A  /  ^  «/<f 


(28) 


This  becomes 


f  (A  -*J\  )  S/*>  &  'd<p'  is.  1 


(29) 


Since,  for  isotropic  scatter,  -** )  is  a  constant 


which  implies 


n+* 

{  cUU<?'  =i  1 


~f  )  —  ,+ji 


(30) 


(31) 


Therefore, 


J-f  <P  (*/  * ')  J  4/1  <f  (*,*•) Uu  ‘ 


(32) 


m  # 

ft 


^yV  *Pi*j 


(33) 


Thus,  the  one-speed,  time  independent,  one-dimensional 
transport  equation  in  a  uniform  and  isotropic  medium  with 


isotropic  scatter  becomes 


(34) 
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The  operator  form  is  the  same  for  all  of  these  simpli¬ 
fied  transport  equations 


£<P-s  =  o 


(35) 


where,  for  the  one-dimensional  case  with  isotropic  scatter, 
the  integro-dif ferential  operator  is 

(36) 


/=  +  Zt  ~  *£■[*' 


Mathematical  Background 

Continuum  problems  often  have  different,  but  equivalent, 
formulations.  In  the  differential  formulation,  the  problem 
is  to  integrate  an  integro-dif ferential  equation  subject  to 
given  boundary  conditions.  In  the  classical  variational 
formulation,  the  problem  is  to  find  an  unknown  function  which 
extremizes  a  functional  subject  to  the  same  given  boundary 
conditions.  A  familiar  example  of  this' dual  formulation  is 
the  equivalence  between  Newtonian  mechanics,  a  differential 
formulation,  and  Hamiltonian  mechanics,  a  variational  formu¬ 
lation  (5:85-105). 

There  are  some  standard  techniques  used  to  translate  a 
differential  formulation  into  a  variational  formulation.  One 
is  to  select  a  quadratic  functional  derived  from  the  differ¬ 
ential  formulation.  The  usual  choice  for  a  quadratic  (6:10) 
is 

(LV,  V)  -  X  (f,  /)  (37) 


where  the  differential  formulation  was 


Lu  -f  =  o 


(38) 


and  with  the  definition  that 


=  [,‘x* 


(39) 


where  signifies  integration  over  the  entire  domain,  and 
y*  is  the  complex  conjugate  of  y.  Since  all  the  variables  in 
the  transport  equation  are  real  valued,  the  complex  conjugate 
notation  will  be  suppressed.  The  integral  to  be  extremized 
then  becomes 


1=  ^[j -lfvfdfj 


(40) 


Another,  but  seldom  used,  quadratic  is 

(L V-fj  Lv-  f  )  (4D 

There  is  one  restriction  on  the  above  approach  —  the 
operator  in  the  resulting  variational  formulation  must  be 
self-adjoint.  Adjointness  is  defined  (8:10)  by  the  property 

[<PX  =  f  t/'fek  (42) 

where  is  the  adjoint  of  the  operator  ^ . 

A  proof  of  this  requirement  follows.  The  transport 
equation  is 


K4>-  S  =  o 


(43) 


and  the  standard  quadratic  is 


Uf,  *)-*■ 

so  the  integral  to  be  extremized  is 

i=  af  <p(x^-2s)^ 

JD 

Taking  the  variation 

SI  =  &j[<PUS-Zs)ei£  =  O 

and 

6l=.Yx\t[s<PiU-Xi)¥ 

using  the  definition  of  adjointness 

SI  =  =  0 

Rearranging  and  combining  terms 

JX-k[/f[W)f-2s}'4  =  0 

or 

but  since  the  variation  is  arbitrary 


•kU*/*)<f-s  =  o 


Clearly,  this  is  equal  to 
is  self-adjoint. 


the  transport 


equation 


(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

only  if 


But 


is  not  self-adjoint  for  the 


transport  equation,  so  something  else  must  be  tried. 


Operator  Formulation 

A  variational  form  of  the  transport  equation  can  be  de¬ 
rived  by  starting  with  a  different  quadratic.  Specifically, 

5  j  (52) 


Now  the  integral  to  be  extremized  becomes 

and 

1=  '/if*  j/V//-  (<*h -s.U?)  +■  s2j 

or 

I 2s&+ 

Using  the  definition  of  adjointness 

X  =  ^jT{ v- s' }& 

Taking  the  variation 

3l*)Ml{4>eZ*-2<PJL*s  t  s*\ek  —  o 


(53) 


(54) 


(55) 


(56) 


(57) 


and 


<fi=  =  o 


(58) 


again  using  the  definition  of  adjointness  and  combining  terms 


yields 


<51=/  WUW-rdik*  o 


(59) 


or 

\£*{£<P-  $)}  dt-o  (60) 

Since  the  variation  is  arbitrary 

/'(jCS-s)  =  ° 

Let 

r=/^-s 

then  there  are  clearly  two  solutions.  Either 

/^-S  -  O  (63) 

which  is  the  solution  to  the  transport  equation  or 

—  o  (64) 

which  is  not. 

Equation  (64)  is  just  the  homogeneous  adjoint  transport 
equation,  which  is  equivalent  to  the  transport  equation  with¬ 
out  sources  and  with  the  velocity  vectors  reversed.  Clearly, 
if  there  is  any  net  absorption,  the  steady  state  solution  is 
identically  zero.  Thus,  this  equation,  applied  to  a  region 
where  there  is  net  absorption,  is  zero  because  is  zero. 


(61) 


(62) 


Therefore,  the  only  solution  to  the  self-adjoint  transport 


equation  is  also  the  solution  to  the  transport  equation. 
As  a  notational  aid,  the  following  definitions  are  made. 


i=r/ 


(65) 


and 

(66) 

Then  the  relation  which  specifies  the  extremization  of  the 
functional  becomes 


Lf-f>=o  <67) 

which  will  henceforth  be  referred  to  as  the  self-adjoint 
transport  equation. 

Integral-Differential  Formulation 

From  the  operator  form  of  the  self-adjoint  transport 
equation,  it  is  a  straightforward  matter  to  derive  the  equa¬ 
tion  in  its  differential  form.  The  operator,  in  the  one¬ 
dimensional  transport  equation  for  the  case  of  isotropic 
scatter  is 


and  its  adjoint  is 


Then,  in  the  self-adjoint  transport  equation 


the  operator 


/y  <?-  ts=o 


L  =  /'/ 


becomes 


+  £-t  L 44  ^  l*1  <^'2 


Collecting  terms  yields 


i=  +•  f_,  <&•' 


(70) 


(71) 


(72) 


(73) 


Thus,  the  one-dimensional,  self-adjoint  transport  equation 
for  isotropic  scatter  becomes 

¥^~[x  [iy (*>*") ‘h”]'*"'  -  <74) 

Since  the  integrals  in  the  double  integral  term  are  indepen- 
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III.  The  Finite  Element  Equations 
Finite  Element  Fundamentals 

The  finite  element  method  is  a  numerical  technique 
which  divides  the  entire  domain  of  the  problem  into  simple 
polygonal  structures  called  elements.  The  solution  over  this 
domain  is  approximated  by  a  set  of  discrete  values,  called 
the  nodal  values,  which  are  located  at  the  nodes,  or  ver¬ 
tices,  of  the  elements.  The  arrangement  of  nodes  and  ele¬ 
ments  is  called  a  mesh.  A  typical  mesh  for  a  two-dimensional 
problem  might  look  something  like  Fig.  1. 


i  €> 


Fig.  1.  Typical  Finite  Element  Mesh 

The  set  of  nodal  values  is  chosen  so  as  to  extremize  the 
variational  statement  of  the  problem.  They  are  computed  by 
solving  the  set  of  simultaneous  equations  derived  from  the 
variational  integral. 

The  solution  is  approximated  within  each  element  by 


interpolating  functions.  If  $  is  the  dependent  variable, 

then  over  the  x-y  plane  <?  is  approximated  by 

* 

<PUj$)  =  %  &L  (78) 

c 

where  n  is  the  number  of  nodes.  This  is  the  finite  element 
approximation.  For  this  element,  <Pc  are  the  nodal  values  at 
the  vertices  of  the  element  and  Mi  are  the  interpolating 
functions.  The  interpolating  functions  are  usually  linear, 
but  higher  order  polynomials  are  sometimes  used. 

In  the  finite  element  approximation,  the  variational 

formulation  has  the  form 

~  k  (79) 

where  b  comes  from  the  source  term.  Because  the  operator 
L  was  self-adjoint,  is  symmetric  and  positive  definite. 

In  the  continuum  formulation,  the  variation  is  taken 
with  respect  to  the  dependent  variable.  In  the  finite  ele¬ 

ment  method,  the  variation  is  taken  with  respect  to  the 
nodal  values.  In  summation  notation,  the  integrrl  to  be 
extremized  is 

I=HS  f  (80) 

1/  it  * 

Taking  the  variation  with  respect  to  the  i'th  node 

<fl=  <8i) 

Since  is  symmetric  and  since  the  indices  are  dummy,  this 


becomes 


Jl  =  Mcj  €'  f  |  ^  ~  ^  ( 82 ) 

thus, 

7^  =  I  *>«  &  ~  k-  =  O  (83) 

J 

or 

(84) 


These  are  the  simultaneous  equations.  The  solution  to  them 
yields  the  nodal  values. 

As  in  any  differential  equation,  the  solution  in  the 
variational  formulation  is  not  unique  until  boundary  condi¬ 
tions  are  specified.  In  the  variational  formulation,  the 
boundary  conditions  are  specified  by  the  requirement  that 
the  variation  be  zero  for  any  points  under  boundary  condi¬ 
tions.  When  the  problem  is  posed  in  finite  elements,  some 
of  the  nodes  are  necessarily  on  the  boundary  and  no  variation 
is  allowed  at  these  points. 

Numbering  the  nodes  such  that  all  the  interior  nodes 
come  first  allows  the  variational  formulation,  in  the  absence 
of  sources,  to  be  restated  as 


i  =  a  ii  1 1] 


£4 1/ 

t 

i 


(85) 


where  s£  is  the  vector  of  interior  nodes  and  is  the  vector 


v.v 
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of  boundary  nodes.  Thus,  the  variation  becomes 


<51  =  Vz  [d?C4„  +  £L,  4?  ft  +  \-°  (86) 

or,  since  the  matrix  is  symmetric, 

JI  ~  &  +~  &i2  H  =■  O  (87) 

and  the  simultaneous  equations  become 

M  (88> 

The  variational  principle  is  good  over  any  sub-domain  of 
the  problem  domain.  Specifically,  it  is  good  over  any  ele¬ 
ment.  Thus, 

<Pl  (89) 

ij 

is  still  valid,  where  £  is  the  number  of  nodes  on  the  ele¬ 
ment  and  is  called  the  local  matrix. 

The  global  matrix  is  created  by  assembling  the  local 
matrices  from  every  finite  element.  The  process  of  assembly 
has  three  steps: 

(1)  Zero  the  global  matrix. 

(2)  For  each  element  in  the  local  matrix,  trans¬ 
late  its  local  indices  into  the  global 
indices,  and  add  it  to  the  proper  global 
element . 

(3)  Do  step  two  for  every  finite  element. 

Since  each  node  is  associated  with  only  a  few  elements,  the 
resulting  global  matrix  is  sparse. 
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One  of  the  strengths  of  the  finite  element  method  is  the 
ability  to  place  elements,  and  therefore  nodes,  where  de¬ 


sired.  In  many  numerical  techniques,  such  as  finite  dif¬ 
ference,  the  nodes  must  be  in  a  rectangular  grid,  and  usually 
with  uniform  spacing.  In  these  other  techniques,  a  lot  of 
computation  must  be  done  in  regions  of  low  interest  just  to 
get  sufficient  resolution  in  the  area  of  high  interest. 

In  the  finite  element  method,  new  elements  can  be  added 
to  the  mesh  arbitrarily.  The  decision  on  where  to  place 
them  is  usually  dependent  upon  how  the  solution  looked  with 
the  original  mesh.  For  instance,  areas  where  the  solution 
is  rapidly  varying  are  likely  candidates  for  a  finer  mesh. 
Since  the  finite  element  method  came  from  a  variational 
formulation,  another  guide  to  placing  new  elements  is  the 
value  of  the  variational  integral  over  each  element.  This 
value  is  also  known  as  the  penalty,  and  the  associated 
variational  integral  is  the  penalty  function.  Elements  with 
the  highest  value  are  sub-divided  into  new  elements  and  the 
problem  is  re-solved.  This  process  can  be  automated  with  a 
computer  routine,  refining  the  mesh  until  some  predetermined 
convergence  criteria  is  reached. 

The  most  commonly  used  shape  for  an  element  in  the  plane 
is  a  triangle.  It  is  the  element  used  in  this  study,  and  a 
brief  review  of  its  properties  follows. 

The  finite  element  approximation  is 

=  £  Mi  <?i 
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,*»  .•*  .v 


(90) 


x  -.w 


The  linear  interpolating  function  is,  over  the  x-u  plane, 


Mi  — 


+  l>;  M  r  a 


Where,  with  i,  j,  and  k  as  cyclic  indices, 


Cic  =  M-  -  Mk 


(91) 


(92) 


be'  =  \  ~ 


(93) 


ti  ~ 


2A-  Xi  X- 

M-  AA:  Mt 


(94) 


(95) 


There  is  a  natural  coordinate  system  for  triangles,  called 
the  triangular  coordinate  system  (3:89).  It  is  depicted  in 
Fig.  2.  The  plane  coordinates  are  related  to  the  triangular 


coordinates  by 


J  l  l  I  Li 

X  _  X,  Xi  Xj  L' 

M  M,  X*)  lL* 

<p,  <  *> 


Note  that  the  relationship 


<P=  <f,h  +  %  U  +  L3 


(96) 


(97) 


implies  that,  for  linear  interpolating  functions, 


Li~M 


(98) 


The  formula  for  integrating  over  the  entire  triangle  is 


f,r  ,1  f  J.  _  ,_r(  j j  rL  _ 
J  L,  Le  L3  eft  -  (  r  +  2). 


(99) 


Fig.  2.  Natural  Triangular  Coordinates 


Two  useful  texts  are  The  Finite  Element  Method  for 
Engineers.  by  Kenneth  H.  Huebner  (7),  and  An  Analysis  of 
the  Finite  Element  Method,  by  George  Fix  and  Gilbert  Strang 
(6).  The  first  text  concentrates  on  applications,  the  second 
one  on  theory. 


The  Variational  Integral 

In  this  study,  the  finite  element  equations  for  triangu¬ 
lar  elements  were  derived  for  the  homogeneous  self-adjoint 
transport  equation  without  sources.  The  integral  to  be  ex- 
tremized  for  this  equation  is 


(100) 


or 


•h  .Si  <S)  f*', . 


/ 1 


»  i 


j 


nnn 


Except  for  the  limits  of  integration,  this  is  also  the  inte¬ 


gral  to  be  extremized  over  each  element.  Specifically, 

z-'a{\ *-  a  <t> iy  <*,«■) « (»' 

+  -  gfr '  J  ^  (102) 

where  M  specifies  integration  over  the  area  of  the  element 
and  is 

dA  “  ^  Am  (103) 

The  next  step  is  to  substitute  the  interpolating  func¬ 
tions  for  (p  .  It  is  desirable  to  use  linear  interpolating 
functions  for  their  computational  efficiency.  However,  the 
class  of  admissible  interpolating  functions  is  limited  by 
the  highest  occurring  derivative  in  the  extremization  inte¬ 
gral  (7:79).  Because  of  the  second  derivative  in  the  first 
term,  linear  interpolating  functions  will  not  converge  to  the 
correct  answer.  An  integration  by  parts  nicely  resolves  this 
difficulty. 

Taking  the  term  with  the  second  order  derivative 

I,  ■  /i  [\ [-*'<!?  Jr  &J]  M  <  104  > 

and  making  the  following  definitions: 

U>=<P  (105) 

thus, 

(106) 
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Therefore, 


@ 


»  f 


x,  =  -%  !*'[<?#  ('  -  j&M  *• 


(107) 


In  a  variational  approach,  the  variation  is  restricted 
to  being  zero  on  the  boundary.  Thus,  the  first  term  in  the 
integration  by  parts,  called  the  surface  term,  is  clearly 
zero  whenever  Y,  and  X*.  are  on  the  boundary.  Since  this 
integral  is  over  a  single  element,  and  since  many  elements 
are  interior  to  the  region,  Y(  and  will  not,  in  general, 
lie  on  the  boundary.  However,  the  value  of  the  surface  term 
has  equal  magnitude  and  opposite  sign  for  adjacent  elements. 
Consequently,  when  the  global  matrix  is  assembled,  all  the 
contributions  from  this  term  cancel.  Thus, 


i=y2  U&fcM 


(108) 


It  is  convenient  to  split  the  extremization  integral 
into  five  integrals,  one  for  each  term.  They  are 


(109) 


which  will  be  referred  to  as  the  streaming  term. 


(110) 


which  will  be  referred  to  as  the  absorbing  term, 


(in) 


which  will  be  referred  to  as  the  first  scatter  term, 


-~vl-  -  ■  'v'l-.'--'  i.'-'-V.t 
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(112) 


which  will  be  referred  to  as  the  second  scatter  term, 


Ts^-Yih' m 


which  will  be  referred  to  as  the  third  scatter  term. 


(113) 


There  is  an  important  distinction  to  be  drawn  among 
the  five  integrals.  The  first  two,  the  streaming  and  absorp¬ 
tion  terms,  are  local,  i.e.,  the  integration  is  over  only  the 
immediate  element.  In  the  scatter  terms,  however,  the  addi- 
tional  integration  over  yOt  implies  involving  other  elements. 
These  terms  are  non-local. 

Most  problems  in  which  the  finite  element  method  is 
used  have  only  local  terms.  The  resulting  global  matrix  is 
sparse  and  only  one  local  matrix  per  element  is  needed.  The 
resulting  computational  efficiency  usually  more  than  compen¬ 
sates  for  the  relative  complexity  of  the  finite  element 
approach. 

In  this  problem,  the  situation  is  not  so  favorable,  but 
since  the  non-localness  is  inherent  in  the  tranport  equation, 
other  numerical  techniques  also  face  similar  difficulties. 
In  general,  this  property  of  the  transport  equation  is  a 
difficulty  in  any  numerical  technique.  How  it  is  handled  in 
this  approach  will  be  described  in  the  section  on  non-local 


■  »  •  *  »  •  •  *_•  •  •  V-  •••••■,  i j,  •  ,  <  ^  ■>. 


terms 


The  Local  Terms 


Starting  with  the  streaming  term 


=  Ua 


(114) 


and  substituting  for  ^  the  finite  element  approximation 


with  linear  interpolating  functions 

3 


< *Cx,  m)  a  $  fiA/ifaM) 


(115) 


The  term  becomes 


(116) 


or 


and 


(117) 


(118) 


Clearly  then, 


(t) 


M 


?  «  Mx^  &  <Ja 


J*  Jx. 


(119) 


where  <**«  is  the  contribution  to  the  local  matrix  from  the 
streaming  term.  The  interpolating  functions  are 


A/t  —  XT  (*lX  ¥  *  tl ) 


Substituting  and  taking  the  derivative  yields 


(120) 


This  integral  can  be  computed  from  the  triangular  integra¬ 
tion  formula.  Substituting  for 


+  A4X  Le  +  Lj 


(122) 


and  integrating  yields 


jjifetA  =  -4~  £#i fa (m,  (U4S  )J  (123; 


+-  M,  fa+M} )  (M3  f^jj  (124) 


thus, 


where 


to  P 
/ yi w  =.  — 1— 


K 


(125) 


(126) 


and  the  i,  j,  and  k  are  cyclic.  Another  form,  which  is  opti¬ 


mized  for  computational  efficiency,  is 


(0  , 


Tva  Viva^ 


(127) 


This  is  the  form  used  in  the  computer  program.  A  copy 
appears  in  Appendix  B. 

Similarly,  where  the  absorption  term 


(128) 


becomes 


J2=  '/A<P;  Z]UiA/;Jj 

t  i  J 


(129) 


but  a  shortcut  can  be  taken  by  recognizing  the  equivalence 
between  the  linear  interpolating  functions  and  the  natural 


triangular  coordinates 


Mi  =.  'jjj-  (acx  +  di)  =  L i 


(130) 


Substituting  the  natural  triangular  coordinates  for  the 


linear  interpolating  functions 


Mij  =  J  Li  Lj  dA 

From  the  triangular  integration  formula 

f  LiLi^M  ~  ~  J  Li  IjdA 


(131) 


(132) 


Thus, 


local  matrix  from  the  absorption  term  is 


.  UJ 

Mi; 


A 

=  ir^t  ' 


I  I  2- 


(133) 


The  Non-Local  Terms 

Before  considering  the  scatter  terms  in  detail,  the 
issue  of  their  non-locality  must  be  dealt  with.  In  Fig.  3, 
a  potential  mesh  over  the  x-u  plane  is  displayed.  The  shaded 
triangle  is  the  current  local  element  and  the  region  between 
the  dotted  lines  is  the  area  of  integration  specified  by 

cfy  eUS. 


'.a. 


The  integration  over  involves  many  elements  other 

than  the  local  one,  a  decided  disadvantage.  Now,  not  only 
do  distant  elements  contribute  to  the  value  associated  with 
the  local  element,  but  often  only  a  portion  of  the  element 
contributes.  To  account  for  this  partial  contribution  would 
seem  to  require  the  existence  of  virtual  elements  (the  parts 
of  the  real  elements).  Their  existences  would  be  temporary 
and  their  structure  would  be  different  for  every  local  ele¬ 
ment.  In  addition,  there  are  now  many  nodes  involved  (in¬ 
stead  of  just  three),  so  the  sparseness  of  the  global  matrix 
is  greatly  reduced. 

In  order  to  avoid  these  difficulties,  a  restriction  on 
the  structure  of  the  mesh  is  imposed.  If  the  numbered  nodes 
from  Fig.  3  were  moved  to  line  up  on  the  dotted  lines  as  in 
Fig.  4,  then  the  integration  over  AA*  involves  only  complete 
non-local  elements.  This  greatly  simplifies  the  derivation 
of  the  finite  element  equations.  In  addition,  the  number  of 
nodes  involved  is  somewhat  less,  so  the  loss  of  sparseness  in 
the  global  matrix  is  less  severe.  In  this  particular  exam¬ 
ple,  nine  additional  nodes  are  involved  in  the  integration 
✓ 

over  M  before  the  numbered  nodes  are  moved,  but  only  seven 
afterwards. 

The  proposed  restriction  on  the  mesh  structure  is  to 
align  the  elements  into  columns  as  in  Fig.  5.  This  will  be 
referred  to  as  the  columnar  scheme.  While  there  is  some  loss 
in  flexibility  of  placement  of  elements,  the  figure  shows 


that  mesh  structures  are  possible  which  still  concentrate  the 
elements  in  chosen  regions. 

✓ 

The  integration  over  the  entire  range  of  XL  is  most 
conveniently  done  one  non-local  element  at  a  time.  The  three 
scatter  terms  become 

—  / 

I  z  j (134) 


Zv  =  <135> 

IV  =  -X  <136) 

—  /  / 

where  XA  is  the  upper  edge  of  the  non-local  element  and  M  is 

—  *  0 

the  lower  edge.  Note  that  both  /A  and  £  are  functions  of  x. 
These  integrals  are  then  evaluated  for  every  non-local  ele¬ 
ment  in  the  column  containing  the  local  element. 

Since  both  the  local  and  non-local  elements  may  point 
in  either  direction,  there  are  clearly  four  cases.  The  cases 
and  the  relevant  coordinate  system  are  summarized  in  Fig.  6. 

The  finite  element  approximation  for  the  three  scatter 
terms  is  a  combination  of  the  forms  used  in  the  streaming  and 
absorption  terms.  For  , 

3 

dp(*,M)  =  5  L-  (137) 

4 

while  for  <p  , 


(138) 


where 


Mi  =  ~J7  (<*/*  +  kj  M'  f  tj  ) 

with 

dj  ~  Mk-  m} 

* 

k  =  X;  -  xk 

Cj  ~~  XcM.k 

and  where  A  is  the  area  of  the  non-local 

P*'_  __L  ' 

*x  “  xd  ^ 

For  the  three  scatter  terms 

X,  =  ->i*  (?(*,*) /V^  M.)  oLu'dti 
2*  =  '/z~-[ <?%  M)J“  M  ifkdil  Jjj'dA 


element 


The  following  finite  element  approximations  are  made 


where 


Thus, 


%  =  jy  (ajx  +  l>J  V+  i}) 


<?x  ZA'  ■> 


(149) 


(150) 


Mcj  =  -  °<  X^7  j  Li  I  ,  (*/*  1-  ^  O  ^  ^ 


Therefore,  the  matrices  for  each  of  the  three  scatter  terms 
become 


(151) 


(152) 


(153) 


M  <  |  f  (* 

/M{;  =  ~  Xf "j  Li]  , 


{'= --t 


Integrating  over 


^  j  +  tj)(*-4')  +  (154) 


MS  = 


2.  ^ 


7  J  Li  UACLj  cM 


(155) 


(156) 


i  "  ^  ZJ\'  J  )  If-J*  +  +  T bj(M  i  sU'J]  <£) 


_  £<  »  f 

X  ZA‘  /  J 


(157) 


(158) 
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(159) 


m 

-Lc 


2 


tof  —  -  fu i {£'+ s’)<=&a 

The  upper  and  lower  limits  of  integration,  fi'  and  44', 
must  be  evaluated,  depending  upon  the  case.  For  Cases  1  and 
3,  the  lower  limit  )  is  specified  by 


thus. 


L2  -  o 


1  1  1  c 

*  *>  }  L*  =  O 

A  Mi  A 


(160) 


(161) 


for  this  to  be  equal  to  zero,  the  determinant  must  be  equal 


to  zero.  Thus, 


I  I  I 

*  b  v  -o 

m'  ui  a' 


(162) 


which  implies 


^  - ^3  (^  -  y,)  +  ^; (y  - ?>)  »  o  (163) 


«  ^  ~L4~-~U  V  +  - - - - - 

f*  S' 


(164) 


However,  x  and  y  have  the  same  range  and  sense,  and  since 


j, '  _  Mi~M'  . 


(165) 


•«  ••  , 


By  symmetry,  the  upper  limit  is 


_  /  M3~M,  -  Mi  y, 

M  ~  ^  +  ?t~>f 

For  Cases  2  and  4,  the  lower  limit  is  specified  by 


i-i  =  0 


thus , 


I  I 


¥»  X*  »  =  <2 


which  implies 


Mb*-}  (?-?*)  =  ° 


/  >4/*  -At i  ,  . 

However,  x  and  y  have  the  same  range  and  sense,  so  this 


just 


M  ~  -  ^  ■+  - - - 


So  again,  by  symmetry,  the  upper  limit  is 

77  '  -  y  1  .flgi  i£k&  (1 

"  >I-y,  >*->/ 

_ / 

In  the  local  matrix  for  the  scatter  terms,  and  ^ 
appear  in  only  two  combinations 


V  -1 


T-  r  •■  >.  tts-ttstt  17 17; 


and 


- '  / 

M  +  41 


(174) 


The  first  combination  becomes,  for  Cases  1  and  3, 

UA.” -id  —  ^  ^  ~  ~^')X  ~  (175! 


or 


^  )] 


(176) 


which  can  be  simplified  through 


M-<U’=  T  1  Jt  [(*3- ?,)](*-*>) 


(177) 


to 


A'-*’  -  --j£(x-y.) 


(178) 


For  Cases  2  and  4, 


it-V=  (4y,  -^%yivr^)]  (179) 


as  before,  this  is 

__  / 


2/J 


^  =  -^r(*-/) 


(180) 


Therefore, 


„  /  / 


=  An  ^(x-y,) 


(181) 


where  is  negative  one  for  Cases  1  and  3,  and  positive  one 
for  Cases  2  and  4. 


41 


2  'V‘2-'2'’ pi  j 


The  second  combination  of  the  integration  limits  is 


O 


JA  -H  UA 

Again,  considering  Cases  1  and  3  first,  this  becomes 

~  )x-h ("/-'*') x  t  -*>?,) i- 

or 

U>-m'  = 

Let 

jt  — *  (-XXj  i  t  -2  aa,) 

and 

k  —  (^j  -f-sM,  —2-m,  ) 

Thus, 

M+-y'=  -  “r  [£x  -  kj//) 

For  Cases  2  and  4, 

*  **'  =  X  t  [4'4)x  +(*yt  ~M#,)+  (Vy, 

or 

Mt-M  -  ^  X  +  (M>  fX4'1 

and  making  the  same  substitution 


(182) 


(183) 


(184) 


(185) 


(186) 


(187) 


(188) 


(189) 


isessssss  mum  &&&&  msesss  ismm  ssssssr 


—  I 

AA  +  Ad 


b<) 


i 

< 

3 

* 

S 

\ 

'« 

(190)  1 


Thus, 

A  +  ~  b'^ 


f 

I 

(191)  i 


Therefore,  the  matrices  for  each  of  the  three  scatter 
terms  become 

m\j  =  -4  ^3  jju  (x-yjfet/x  +  +2-^1;  (Ax  -  ky,)']  M  ( 192 ) 

Mcj  =■  'jzfliix-J/i)  <2jM  *64  (193) 

Alcj  =  ”  ^  j L^i  T  (Ax  -  ky,)  cl  A  (194) 

Note  that  there  is  a  lot  of  parallelism  among  these 
three  expressions.  Later  mathematical  manipulations  can  be 
reduced  by  a  judicious  rearranging  of  the  terms  within  each 
expression.  Thus, 

A>ij  =  -  <K 71, > 

-^7-  Kf  u[£**~  (H  +  h,)  7  *  ^  <195> 

and 

=  "f1'  ^  IT  Q fh(x-%)AAd4  (196) 


and 
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(197) 


a 


/.  ** 
vv 


O 


* 


M 


u 


=:  -  -|*  2^1  fk  \j-i7~  UikJfyx  f  kyf'Jo/A 


Before  proceeding  to  evaluate  each  scatter  term,  it  is 
convenient  to  build  a  table  of  integrals.  The  five  inte¬ 
grals,  which  are  based  upon  the  triangular  integration  formu¬ 
la,  are  derived  in  Appendix  A  and  are  listed  in  Table  I. 

A  note  about  the  notation  used  in  Table  I  is  in  order. 
The  vector  of  interpolating  coefficients  is 


L 


Li 

Lj 


=M 


(198) 


The  area,  A  ,  is  for  the  local  element. 

Because  of  the  parallelism  in  the  form  of  the  scatter 
terms,  it  seems  most  convenient  to  evaluate  the  simplest 
first.  Starting,  then,  with  the  second  scatter  term,  its 
matrix  can  be  rewritten  as 

W 


M 


—  2*  -  ft  J 


(199) 


which  becomes 


m g  =sr%>7k> 


(200) 


For  the  symmetric  cases,  Cases  1  and  4,  where  s=  X, 

<  _  A 


~ 2*  bod1  ^  f— ' " ^a) -M. 


(201) 


or 


44 


(202) 


( 


K‘  =  ^  /?'!  aJf< &-*■)&>} M 


Thus, 


lH>:  =  *•/?,,  -JZ'  %%><(■&!,  M- 


2.  /c"  (,0d 


(203) 


which  is 


,w  _  J&L-ii  a’. 


Mi  - 


<0  ~  2.  (JA 


a,mkM 


(204) 


Explicitly , 


*  s  A  I*  ’  H  ' 

a/=- f-^7  3  f  *  -»•  * 

L3  v  fJhJ 

Similarly,  for  the  antisymmetric  cases,  where  y,  s  X, 

<+C  =  ^£7>a  ZJ1  %'{(*,-*>)&<$  4 


(205) 


K;~  it  *11  ZZi  aJ  *>•  4L 


’w  3.  LOi 


Thus, 


j  .  W  $*5  ^  / 

^  Z2  4  »-* 


Explicitly, 


.  CV 
Mr; 


Sx  jL 


6  2  2 
2  2  1 
2  /  2 


*  K* 


(206) 


(207) 


(208) 


(209) 


£3 


«W»«5®»3 
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Evaluating  the  third  scatter  term  next,  its  matrix  can 
be  rewritten  as 

=■  -  IT  d]  fir  L  (£<z- %x  +  ky,1)  JA  (210) 

where 

£  ~  +  k)  (211) 

Expanding 

A  '(%  +  ?*)  (212) 

Substituting  for  and  ,  and  evaluating  the  integrals 

yields 

^6'  4}T  (-2a' [(*,&, M bsu)*- fy  yy&j] 

[lx, -  Xy,  +y,*&*x] 

(213) 

/ 

The  coefficient  of  M,  is  independent  of  case.  Working  with 
just  that  term,  it  becomes,  for  all  cases 

«2V[-  Ugt+^toi)x+- (x,  +x>)m* x  - tei &+JQ 

Rearranging  and  combining  terms 

2*!  [*,&>£  f 


(214) 


(215) 


Expanding  the  matrices,  and  then  contracting  them  on  the 
basis  of  X9-Xj  yields 


3  n 

3  12 


<*  * 
X  3 
Z  3 


Rearranging  around  the  quadratic  terms  in  x  yields 


f  14  x,?~  % 

x  <  3  y,2  —  £  y,  yA  4-  3  x  * 
/  3V  -  £  y^2  ^  ~hXt 


(217) 


But  this  is  just 


u:  3  ^ 


(218) 


The  coefficients  of  Mt  and  M3  are  identical.  There¬ 
fore,  only  one  needs  to  be  evaluated.  However,  they  are  case 
dependent.  Thus,  for  the  symmetric  cases 

I~]  (219) 


+  x,w*(x,z-x)] 


(220) 


Expanding  the  matrices  as  before 


I  *  (*  I  *>  I  tO  10  I*,- *» 

3  12  W  +  y,  5"  ir  [*,  -  *i 
I  3  t*  S  IS 


(221) 


Rearranging  around 


vv 


[H*,-*-  Lxt  -10/,)  ^  (  Ui-x<)c  (* 

£x2--*,)C?x,  +  /2Xj -/**.)  r  —  {Ur*'fu 

-us^-w,)  J 


(*>-*,)  >2 


which  is 


d 


z 


Similarly,  for  the  antisymmetric  cases 


or 


j\x,  ’f’  ^■0 


Again,  expanding  the  matrices 


U>  •<£*,-*) 


L  v 
z  z 

Z .  3 


[; 


+  X2 


)o  hT\ 
S  15 
S’  IS- 


X,-*f 

Xj  ~  Xz 


Rearranging  around  X,  —  X2  yields 


which  is  just 


'(X,  -X*)  (£x,  +  ¥Xj  -/OX«) 

(X,-X»)  (ZX,  +  ?X2-  5X,J 
(X,  ~  X*J  (2X,  +  3X2  -  5Xd , 


.2^/ 


3 

L'J 


dx 
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Therefore,  for  the  symmetric  case 


.  ^ 

= 


Si.  A, 

2  (Jet 


if  3  3 


3  (*  (*  JA 

ILL 


(?)  <,  A  / 


where 


&,=  3  t  * 

j  U  t» 


V  3  3 


And  likewise,  for  the  antisymmetric  case 


/%•  = 


<  A  A1*  3  3 

km*  *  >  ’  J 


3  /  / 


»?«  -  St.  * 
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where 


if  3  3 

&=  3  J  / 

3  I  I 


(229) 


(230) 


(231) 


(232) 


(233) 


(234) 


The  first  scatter  term  is  most  easily  evaluated  one 


piece  at  a  time.  The  first  piece 


(235) 


vV.V.V.V.V 
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—  <5(  7?/? 
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LOdx 


For  the  symmetric  cases 


lO<C 


10  5  F 

y  io  5 

5  S  \o 


~  */ 

*5  -  X, 


or 


or 


A  J  1 

”<?  y  f~ 

1 

ro 
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1 

• 
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u 
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For  the  antisymmetric  cases 


*  ^  m'  tj 


y 

y” 

y 
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5 

f 

L 

F 

/* 

or 


<*n  u<e  LJ 


)D  0  o 
y  o  o 
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hyJ 


or 

A 

cod 


The  second  piece 


(243) 


by  symmetry  with  the  second  scatter  term  is 

A 


■~(k 


LOct 


C  3  3] 

V 

?  9  ¥ 

3  H  t  \ 
-  J 

U_ 

a\ 


or 


-*Tu 


for  the  symmetric  cases,  and 

-"S? 


z  z~ 

\ 

2  1  1 

2  I  2 

A 

4: 


or 


for  the  antisymmetric  cases. 
The  third  piece 


I 


(244) 


(245) 


(246) 


(247) 


(248) 


by  symmetry  with  the  third  scatter  term  is 


(250) 


A  \ /  * 


for  the  symmetric  cases,  and 


^  i  *  3  3  I 

*  ted  V  3  1  i  4 

3  I  ' 


(251) 


,4  .  ' 

-4  ^  V  ®j*t 


(252) 


for  the  antisymmetric  cases. 


Putting  the  three  pieces  together  for  the  symmetric 


cases  yields 


$  =  -+&e4l>s  +  h>Itj] 


(253) 


Expanding  the  matrices 


(*  ,  a  ***■<] .  ^ ,  ]°  *  flfo ) 

^’""'^55  1  J  ?  ^  HU  *s  If +  0/0  S’  I  tjk  254) 

3  V  ?  W  3  4  U  U'l  J  [0  r  /q[»J  / 


Since  =  -Vj  ,  the  lower  right  four  elements  in  the  first  ma¬ 


trix  can  be  rearranged.  Similarly,  all  of  the  elements  in 


the  last  matrix  can  be  freely  rearranged.  Thus, 


4  1^3  3 


,[4  3  3~)M  ,  r*3  3")fl 


3  6  6  *>  »l  (  } 


SFI 


(256) 


The  bracketed  term  can  be  evaluated  by  recognizing  that  the 
natural  triangular  coordinates  in  the  non-local  element  are 


l~j  =  kj m'+  t-'] 


(257) 


which,  for  a  point  at  its  corresponding  node,  is  equal  to  one 


I  -  7~7  t  Cj'J 


(258) 


but  at  any  other  node  is  equal  to  zero 

<5>= 

Thus,  the  bracketed  term  is 


iV,j 


(259) 


(260) 


and  therefore,  for  the  symmetric  case,  the  matrix  for  the 
first  scatter  term  is 

A  A' 


(261) 


Putting  the  three  pieces  together  for  the  antisymmetric  cases 
yields 

O)  A 


ev^'  */*■  h>  i  ti  ] 


(262) 


Expanding  the  matrices 
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zr 


)0  O  O 
TOO 


TOO 


t;  y  (263) 


While  it  is  still  possible  to  rearrange  the  terms  in  the 
matrices  corresponding  to  5,  it  does  no  good  because  the 

term  is  no  longer  equal  to  the  -M,  term  in  the  second 
matrix  as  it  was  in  the  symmetric  cases.  Additionally,  this 
author  was  unable  to  simplify  this  expression  to  any  signifi¬ 
cant  degree. 

A  simple  answer  was  obtained  by  a  different  derivation 
and  it  is  included  in  Appendix  G.  However,  it  does  not  build 
on  the  derivation  presented  here.  Nonetheless,  it  does  ob¬ 
tain  a  useful  answer,  one  which  was  confirmed  independently 
by  Dr.  Shankland,  and  which  is  in  agreement  with  the  answer 
one  might  anticipate  from  the  results  in  the  symmetric  case. 
Assuming,  for  the  moment,  that  this  answer  is  correct,  the 
matrix  for  the  first  scatter  term  for  the  antisymmetric  case 
is 


•—  ^  IQ  U  —  ^ 


(264) 


where 


3  3 
I  I 
/  I 


(265) 
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IV.  The  Numerical  Solution 


Benchmark  Problem 

A  picture  of  the  benchmark  problem  used  in  this  study  is 
shown  in  Fig.  7. 


Fig.  7.  Benchmark  Problem 


The  benchmark  problem  consisted  of  a  mono-energetic 
source  of  particles  which  impinged  from  the  left  upon  a  blan¬ 
ket  three  mean  free  paths  thick.  The  source  was  Lambertian, 
and  constituted  the  first  of  two  boundary  conditions.  The 
blanket  was  surrounded  on  both  sides  by  vacuum  and  an  albedo 
of  zero  was  assumed  so  the  returning  flux  from  the  right  was 
also  zero.  This  was  the  second  boundary  condition.  The 
total  macroscopic  cross-section  was  chosen  as  unity  so  that 
it  could  serve  as  a  scaling  factor.  The  scattering  macro¬ 
scopic  cross-section  was  either  zero,  which  eliminated  the 
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scattering  terms,  or  half  the  total  cross-section.  Both 
macroscopic  cross-sections  were  constant  throughout  the 
blanket.  There  were  no  inbedded  sources  in  the  blanket. 


Algorithm 

The  algorithm  for  the  benchmark  problem  follows  the 
usual  finite  element  approach,  except  for  the  handling  of  the 


scatter  terms.  Specifically: 


(1)  Partition  the  global  matrix  into  four  sub¬ 
matrices  according  to  the  interior  and 
boundary  nodes. 

(2)  Zero  the  global  matrix. 

(3)  Define  the  mesh. 


(4)  Compute  the  contribution  to  the  local  matrix 
from  the  streaming  term  and  assemble  it  into 
the  global  matrix. 

(5)  Compute  the  contribution  to  the  local  matrix 
from  the  absorption  term  and  assemble  it  into 
the  global  matrix. 

(6)  Save  the  contributions  from  the  streaming  and 
absorption  terms. 


(7) 


For 

(a) 


(b) 


(c) 


(d) 


each  non-local  element  in  the  column: 
Compute  the  contribution  to  the  local 


matrix  from 
assemble  it 
Compute  the 
matrix  from 
assemble  it 
Compute  the 
matrix  from 
assemble  it 
Save  the  contributions 
scatter  terms. 


the  first  scatter  term  and 
into  the  global  matrix, 
contribution  to  the  local 
the  second  scatter  term  and 
into  the  global  matrix, 
contribution  to  the  local 
the  third  scatter  term  and 
into  the  global  matrix. 

from  the  three 


(8)  Compute  the  force  vector. 

(9)  Solve  the  system  of  simultaneous  equations. 
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(10)  Compute  the  element  penalties  and  the  global 
penalty. 

(11)  Print  the  results. 

Implementation  Notes 

The  computer  program  was  written  in  FORTRAN  77,  under 
the  UNIX  operating  system,  release  number  4.1  by  Berkeley, 
and  should  meet  the  77  standard  (2)  in  all  specifics.  While 
it  was  not  optimized  for  computational  efficiency,  it  does 
strive  for  clarity,  and  its  structure  should  aid  attempts  at 
modification  and  improvement.  All  variables  are  passed  in 
the  call  statements;  there  are  no  COMMON  blocks.  All  of  the 
arrays  which  are  dependent  upon  the  mesh  are  dimensioned  via 
parameters.  This  minimizes  storage  space,  but  still  allows 
for  easy  change  in  the  size  of  the  arrays  when  a  new  mesh  is 
used. 

While  comments  are  used  liberally  throughout  the  code, 
the  concentrated  documentation  is  in  the  subroutine  DOCUM. 
It  consists  of  three  glossaries:  a  glossary  of  routines, 
one  of  parameters,  and  one  of  variables. 

The  main  program  is  routine  SATE.  The  first  routine  it 
calls  is  INITAL.  Subroutine  INITAL  creates  the  mesh  by  read¬ 
ing  in  its  definition  from  a  file,  computes  the  areas  of  the 
elements,  and  zeros  the  global  matrix. 

Next,  the  global  matrix  is  assembled.  It  is  stored  in 


four, 

two-dimensional 

array 

variables:  Mil, 

Ml  2, 

M21,  and 

M22. 

The  contributions 

to  the 

local  matrices 

are 

computed 

.  ■  -T  .  -  .  •  .V  _•  .V  ,-W  N  .  .V 
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by  five  subroutines,  one  for  each  of  the  terms  in  the  extrem- 
ization  integral.  The  five  subroutines  are  STREAM,  ABSORB, 
SCAT1,  SCAT2,  and  SCAT3.  The  routine  CASEDT  determines  the 
case  of  the  orientation  of  the  local  and  non-local  elements, 
for  use  by  the  scatter  routines.  The  contributions  to  the 
local  matrices  are  combined  into  temporary  matrices  and  saved 
in  MAT(LOC,I,J)  for  later  use  in  the  calculation  of  the 
penalties.  The  contributions  from  the  streaming  and  absorp¬ 
tion  terms  are  combined  into  one  temporary  matrix.  The  con¬ 
tributions  from  all  three  scatter  terms  are  also  combined 
into  temporary  matrices  but  there  is  now  one  temporary  matrix 
for  each  non-local  element  in  the  column.  This  is  a  memory 
intensive  approach,  but  it  makes  debugging  easier  and  it  does 
avoid  having  to  recompute  these  matrices. 

When  the  global  matrix  is  completely  assembled,  the 
force  vector  is  computed.  It  is  just  the  negative  product  of 
the  M12  matrix  and  the  vector  of  boundary  nodes.  The  Mil 
part  of  the  global  matrix  is  copied  into  the  double  precision 
array  MM  for  two  reasons.  First,  the  implementation  of  IMSL 
on  the  VAX  11/780  is  in  double  precision  only.  Second,  the 
matrix  sent  to  the  IMSL  routine  suffers  a  Gauss-Siedel  decom¬ 
position,  and  is  no  longer  available  for  use  in  computing  the 
global  penalty. 

The  solution  is  computed  by  the  IMSL  routine  LEQIF, 
which  is  a  virtual  memory  version  of  the  linear  equation 
solving  routine. 
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The  penalty  is  computed  in  two  ways.  First,  the  element 
penalties  are  computed.  They  are  useful  as  guides  in  refin¬ 
ing  the  mesh.  The  global  penalty  is  also  computed,  but  for 
comparison  purposes  only,  since  the  sum  of  the  element  penal¬ 
ties  should  be  equal  to  the  global  penalty. 

The  results  are  printed  on  the  standard  output  by  the 
routine  OUTPUT.  It  offers  several  options  for  easy  tailoring 
of  the  desired  information. 

The  last  routine  is  DEBUG.  It  is  a  low  overhead  utility 
routine  which  can  be  used  to  quickly  get  diagnostics  about 
the  program  without  cluttering  it  up  with  print  statements. 
It  is  no  substitute  for  a  symbolic  debugger,  but  it  is  supe¬ 
rior  for  problems  with  a  lot  of  intermediate  data. 

The  major  variables  include  the  previously  mentioned 
arrays,  the  global  matrix  (Mil,  M12,  M21,  M22),  and  the 
record  of  the  contributions  to  the  local  matrices.  In  addi¬ 
tion,  there  are  three  arrays  which  define  the  mesh  in  use. 
They  are  CORDND ,  PTNODE ,  and  AREAS.  CORDND  (coordinates  of 
the  nodes)  contains  the  xm  -coordinates  of  all  the  nodes. 
PTNODE  (pointer  to  the  nodes)  specifies  which  nodes  are  a 
part  of  each  element  by  pointing  to  CORDND.  AREAS  is  just 
the  area  of  each  element. 


V.  Results 


Local  Terms 

The  first  use  of  the  finite  element  equations  was  to 
solve  the  benchmark  problem  for  a  blanket  which  consisted 
of  a  purely  absorbing  medium.  This  tested  the  streaming  and 
absorption  terms  without  the  added  complexity  of  the  scatter¬ 
ing  terms.  The  transport  equation,  with  absorption  only,  is 

+*  —  &  (266) 
This  problem,  when  solved  analytically,  is 
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Evaluating  the  integration  constant,  C,  from  the  boundary 
condition  of  a  Lambertian  source,  yields 

-  4r  X 

-  Ae  M>  O  (268) 

Five  computer  runs  were  done,  one  run  for  each  of  the 
five  meshes.  The  detailed  results  of  these  runs  are  included 
in  Appendix  D.  Some  of  the  notable  features  are  discussed 
below. 

The  analytic  solution,  for  the  case  of  M.  =  1,  is  shown 
in  Table  II,  along  with  the  nodal  values  and  the  global 
penalties  from  the  five  meshes. 

It  can  be  seen  from  Table  II  that,  as  the  number  of 


elements  in  the  mesh  increases,  the  answer  tends  to  converge 
to  the  analytic  solution.  The  difference  between  Mesh  2  and 
Mesh  3  is  not  an  exception,  since  these  two  meshes  have  the 
same  number  of  elements.  The  only  difference  between  them  is 
the  orientation  of  two  adjacent  elements.  (Refer  to  the  mesh 
maps  in  Appendix  C.)  This  small  change  in  going  from  Mesh  2 
to  Mesh  3  produced  a  noticeable  change  in  both  the  answer  and 
the  penalty.  The  answer  is  worse  in  Mesh  3  than  in  Mesh  2 
and  the  penalty  values  reflect  it. 


Table  II 

Nodal  Values  from  Local  Terms  for  -  1.0 


Solution 

Range 

Penalty 

Number 

of 

Elements 

x=0.00 

x=. 375 

x=. 750 

x=l . 50 

x=3.00 

Analytic 

1.00 

.687 

.472 

.223 

.050 

— 

— 

Mesh  1 

1.00 

— 

.316 

.116 

.042 

.300 

9 

Mesh  2 

1.00 

— 

.413 

.157 

.056 

.271 

13 

Mesh  3 

1.00 

— 

.389 

.159 

.057 

.278 

13 

Mesh  4 

1.00 

.637 

.385 

.144 

.052 

.264 

20 

Mesh  5 

1.00 

.657 

.423 

.131 

.028 

.247 

32 

The  inclusion  of  Mesh  3  was  to  test  the  sensitivity  of 


the  answer  to  the  orientation  of  the  elements.  Additionally, 
each  mesh  had  one  element  with  zero  area  as  a  test  of  the 
routine's  ability  to  handle  extreme  cases. 


For  Mesh  5,  for  the  nodes  near  the  region  where  the 
elements  are  concentrated,  the  finite  element  solution  is 
within  10  percent  of  the  analytic  solution.  For  a  mesh  of 
only  32  elements,  this  is  reasonable  accuracy. 

In  Mesh  5,  the  nodal  values  at  x  =  1.50  and  x  =  3.00  are 
both  further  from  the  analytic  solution  than  for  most  of  the 
earlier  meshes.  One  possible  reason  for  this  is  that  the 
style  of  Mesh  5  is  very  different  from  the  other  meshes.  In 
Mesh  5,  the  elements  are  concentrated  in  the  upper  right  cor¬ 
ner.  This  left  some  very  large  elements  in  the  opposite  cor¬ 
ner.  Another  factor  is  that  the  shape  of  the  elements  near 
these  nodes  is  elongated.  For  the  finite  element  method,  the 
preferred  shape  is  a  near  equilateral,  since  this  shape  is 
better  at  approximating  a  solution  with  different  curvatures 
in  each  direction. 

In  Table  III,  the  Mil  portion  of  the  global  matrix  from 
Mesh  4  is  shown.  This  is  the  part  of  the  global  matrix  which 
couples  the  interior  nodes  to  themselves.  The  symmetry, 
bandedness,  and  sparseness  are  clearly  evident.  A  little 
inspection  reveals  that  it  is  also  positive  definite.  These 
are  all  the  features  expected.  The  banded  structure  is  typi¬ 
cal  whenever  an  orderly  numbering  of  the  nodes  is  done  and, 
like  symmetry,  it  can  be  used  to  speed  program  execution  or 
to  reduce  memory  requirements.  However,  in  a  program  which 
uses  automatic  mesh  refinement,  the  bandedness  will  be  dis¬ 
rupted  with  the  addition  of  each  new  element.  In  order  to 


maintain  the  banded  structure,  it  is  necessary  to  utilize  a 
re-numbering  routine.  A  well-constructed  re-numbering  rou¬ 
tine  will  not  only  maintain  the  banded  structure,  but  it  will 
also  minimize  its  "bandwidth". 

Table  III 

Local  Term's  Mil  Global  Matrix  from  Mesh  4 
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.00 

.00 

.03  -. 

00 

.99 

.02 

• 

o 

o 
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00 
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o 
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• 

03 
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o 

• 

• 

o 

to 

.63  . 

65 

.00 

.00 

.03  1. 

00  .00  1.01 


05  1.16 


From  Appendix  D,  it  is  clear  that  in  every  mesh  there 
are  nodes  whose  fluxes  are  negative.  The  magnitudes  are 
small,  but  the  values  are  clearly  non-physical.  The  nodes 
where  these  negative  values  occur  are  expected,  from  the 
analytic  solution,  to  be  small  or  zero,  so  the  difference  is 
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not  great.  Furthermore,  small  differences  are  not  unex¬ 
pected  in  a  numerical  technique,  but  in  some  applications  a 
negative  flux  can  be  very  inconvenient. 


One  potential  solution  is  to  use  a  technique  called 
quadratic  programing  (10).  In  this  technique,  constraints 
are  imposed  upon  the  values  of  the  nodes  so  they  cannot  go 
negative.  Clearly,  this  will  not  improve  the  approximation 
to  the  extremization  integral;  indeed  it  can  only  worsen  it, 
but  the  effect  should  be  small,  and  it  may  be  more  tolerable 
than  having  negative  values. 

Non-Local  Terms 

The  second  use  of  the  finite  element  equations  was  to 
solve  the  benchmark  problem  where  the  blanket  had  a  non-zero 
scattering  cross-section.  The  response  of  the  scatter  terms 
was  the  main  interest. 

Five  computer  runs  were  done,  one  run  for  each  of  the 
five  meshes.  The  detailed  results  of  these  runs  are  included 
in  Appendix  E.  Some  of  the  notable  features  are  discussed 
below. 

It  can  be  seen  from  Table  IV  that,  as  the  number  of  ele¬ 
ments  in  the  mesh  increases,  the  global  penalty  increases. 
These  results  exactly  parallel  the  results  for  the  local 
terms,  right  down  to  the  inversion  of  Meshes  2  and  3.  What 
really  stands  out,  however,  is  the  difference  in  the  magni¬ 
tudes  of  the  penalties  between  the  local  runs  and  the  scat¬ 
tering  runs.  In  the  scatter  runs,  the  penalties  are  lower  by 


a  factor  of  nearly  three.  The  decrease  is  anticipated  since 
the  inclusion  of  the  scatter  terms  smoothes  the  distribution 
of  the  flux.  A  function  with  a  smooth  distribution  is  better 
approximated  by  finite  elements  than  a  function  whose  distri¬ 
bution  is  rapidly  changing. 

Table  IV 

Penalty  Values  for  Scattering  Benchmark  Problem 


Solut ion 

Penalty 

Number 

of 

Elements 

Mesh  1 

.125 

9 

Mesh  2 

.109 

13 

Mesh  3 

.114 

13 

Mesh  4 

.102 

20 

Mesh  5 

.097 

32 

The  penalty  values  for  the  individual  elements  followed 
about  the  same  patterns  in  the  scattering  solution  as  in  the 
absorption  only  solution.  In  particular,  the  relative  fluc¬ 
tuations  from  element  to  element  were  about  the  same  in  both 
solutions. 

The  scattering  benchmark  was  solved  independently  by 
another  numerical  technique  called  Ln  (9).  This  technique 
was  recently  developed  at  the  Air  Force  Institute  of  Technol¬ 
ogy  by  LCDR  Kirk  Mathews,  U.  S.  Navy,  who  graciously  ran  a 
solution  of  the  scattering  benchmark  problem.  The  Ln  results 


served  as  a  baseline  solution  whose  accuracy  was  assumed  to 
be  sufficient  for  the  objectives  of  this  study.  The  output 
from  the  Ln  program  is  the  scalar  flux,  not  the  phase  space 
scalar  flux  that  the  routine  SATE  outputs,  so  some  intermedi¬ 
ate  calculations  were  necessary  before  a  comparison  was  pos¬ 
sible. 

The  intermediate  calculations  were  done  by  hand  and  were 
simple  linear  interpolations  between  the  grid  points  of  the 
Ln  scheme  and  between  the  nodes  in  the  finite  element  scheme. 
A  comparison  of  the  interpolations  is  shown  in  Table  V.  The 
comparison  is  for  Mesh  5  only. 

Table  V 

Scalar  Flux  for  Scattering  Benchmark  Problem 


Number 

Solution 

Range 

Penalty  of 

_  Elements 

x=0.00  x=.375  x=.750  x=1.50  x=3.00 
Ln  .296  .191  .129  .064  .016 

Mesh  5  .285  .235  .182  .091  .036 


The  finite  element  solution  disagrees  with  the  Ln  solu¬ 
tion  by  about  35  percent.  At  x  =  0,  the  value  of  the  flux 
is  driven  by  the  boundary  conditions;  therefore,  the  dif¬ 
ference  at  this  point  is  a  rough  measure  of  the  accuracy  cf 
the  intermediate  calculations.  The  apparent  uncertainty  in 
the  intermediate  calculations  is  clearly  not  enough  to 


explain  the  larger  differences  over  the  mesh  as  a  whole. 

The  finite  element  solution  is  larger  at  every  point, 
which  indicates  a  systematic  error.  It  seems  that  the  scat¬ 
ter  terms  are  over-active;  that  is,  their  contribution  to  the 
solution  is  disproportionate  to  the  contribution  time  from 
the  local  terms. 

Before  addressing  possible  causes  of  the  difference,  the 
global  matrix  should  be  inspected.  The  Mil  portion  of  the 
global  matrix  from  Mesh  4  is  shown  in  Table  VI.  Mesh  4  is 
shown  instead  of  Mesh  5  for  several  reasons.  First  and  fore¬ 
most,  Mesh  4  possesses  all  of  the  significant  features  from 
any  of  the  global  matrices  from  the  scattering  runs.  Addi¬ 
tionally,  the  global  matrix  shown  in  Table  III  for  the  local 
terms  was  also  taken  from  Mesh  4,  and  so  a  direct  comparison 
is  possible.  Lastly,  the  global  matrix  from  Mesh  4  is  the 
largest  from  the  set  of  five  meshes  that  will  conveniently 
fit  on  a  page. 

The  scattering  global  matrix  displays  a  typical  banded 
structure  and  it  is  positive  definite,  both  expected  fea¬ 
tures.  Also,  it  is  much  less  sparse  than  the  global  matrix 
from  the  absorption  only  runs.  This  is  because  of  the  non¬ 
localness  of  the  scattering  terms.  The  decrease  in  the 
sparseness  is  also  an  expected  feature.  The  critical  fea¬ 
ture,  however,  is  the  lack  of  symmetry.  With  the  two  digits 
of  precision  used  in  outputing  the  matrix,  it  is  clear  that 
many  of  the  cross-diagonal  terms  are  not  equal.  While  the 
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differences  are  small,  this  lack  of  symmetry  is  inconsistent 
with  a  self-adjoint  equation.  Clearly,  the  local  matrices 
for  the  scatter  terms  are  being  computed  incorrectly. 

Table  VI 

Scatter  Term's  Mil  Global  Matrix  from  Mesh  4 


hi 

mn 


00  -.65 

.00 

.00  .00 

02  -.00 

-.47 

.00  .00 

03  -.03 

-.18 

-.18  .01 

57  .01 

-.01 

-.24  .00 

00  1.39 

-.03 

-.01  -.65 

01  -.03 

.98 

.07  -.02 

24  -.01 

.08 

.69  -.03 

00  -.66 

-.02 

-.03  1.10 

00  .00 

-.31 

.01  -.05 

00  .00 

00  .00 

.00 

.00 

-.20  .01 

-.00  -.31 

00  .00 

15  .01 

’6  -.06 
>6  .55 

18  -.03 
.2  -.01 


There  are  at  least  three  possible  sources  of  the  error. 
First,  the  local  matrices  may  not  have  been  derived  correct¬ 
ly,  second,  the  finite  element  equations  may  not  have  been 
correctly  coded  into  the  computer  routines,  and  third,  the 
calculation  of  the  self-adjoint  operator  may  be  incorrect. 
If  either  of  the  first  two  possibilities  is  the  source  of 


the  error,  then  correcting  them  should  restore  symmetry  to 
the  global  matrix. 


The  local  matrices  from  the  second  and  third  scatter 

terms  are  not  symmetric  (see  the  debug  output  in  Appendix  F). 

This  was  apparent  at  the  time  of  the  derivation  and  was  not 

of  concern  for  two  reasons.  First,  the  self-adjoint  property 

the  usual  terms  found  in  self-adjoint  equations.  Recall  that 

the  scatter  terms  in  the  extremization  integral  are  to  be 

/ 

integrated  over  the  entire  range  of  M  .  In  addition,  in  this 
derivation  they  are  evaluated  over  one  non-local  element  at 
a  time.  For  these  reasons,  the  contributions  to  the  local 
matrices  from  the  scatter  terms  are  not  expected  to  be  sym¬ 
metric. 

What  should  happen,  however,  is  that  the  combined  con¬ 
tributions  to  the  local  matrix  from  all  the  scatter  terms  be 
symmetric.  Obviously,  this  did  not  happen  since  the  global 
matrix  was  not  symmetric. 

The  possibility  remains  of  coding  errors  in  the  computer 
program.  Several  errors  were  found  in  the  final  days  of  this 
study  and  there  are  no  guarantees  that  all  the  bugs  were 
eradicated. 

In  the  results  for  the  local  terms,  the  fluxes  near  the 
right  boundary  did  not  converge  consistently  as  the  mesh  was 
refined.  It  was  attributed  to  the  relative  coarseness  of  the 
meshes  in  that  region.  It  is  possible,  however,  that  another 
effect  is  in  operation. 
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When  the  self-adjoint  operator  was  derived,  it  was  com¬ 
puted  from  ft. .  If  instead,  it  is  derived  from  the  ex- 
tremization  integral  (Equation  53),  then  a  surface  term 
appears.  The  surface  term  does  not  change  the  equations  of 
motion,  but  it  does  prescribe  that  the  conventional  trans¬ 
port  equation  be  satified  on  the  free  boundaries,  a  proper 
result.  This  term  is  missing  in  the  derivation  done  in  Chap¬ 
ter  II  and,  therefore,  there  may  be  an  incorrect  boundary 
condition  in  effect.  This  would  put  strains  on  the  solution 
and  would  seriously  perturb  it  near  the  free  boundaries.  It 
would  be  useful  to  have  a  uniform  mesh  in  investigating  this 
possibility. 

In  the  scattering  runs,  not  a  single  node  was  found  to 
have  negative  values,  not  for  any  of  the  meshes.  This  is 
reasonable,  since  the  process  of  scatter  tends  to  smooth  and 
spread  the  distribution  of  the  flux. 

The  real-time  run  time  for  the  benchmark  problem  with 
scatter,  using  Mesh  5,  was  7  seconds.  This  was  on  a  VAX  11/ 
780  with  no  other  major  programs  running.  The  FORTRAN  com¬ 
piler  on  this  particular  machine  is  very  inefficient,  how¬ 
ever,  since  the  FORTRAN  is  translated  into  the  programing 
language  C  before  compilation.  In  addition,  this  implementa¬ 
tion  of  the  finite  element  equations  was  not  optimized  for 
computational  efficiency. 
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VI.  Conclusion  and  Recommendations 

Conclusion 

This  study  accomplished  three  major  tasks  toward  a  gen¬ 
eral  solution  of  the  transport  equation  using  finite  ele¬ 
ments:  the  recasting  of  the  transport  equation  into  a  self- 
adjoint  form,  the  derivation  of  the  finite  element  equations 

<r. 

for  the  one-dimensional  case  with  isotropic  scatter,  and  the 
implementation  of  these  equations  in  a  computer  code.  While 
the  solution  has  the  expected  features,  some  questions  remain 
about  residual  errors  in  the  scatter  terms. 

Several  features  of  the  finite  element  solution  stand 
out.  The  integer  nature  of  the  cores  of  the  local  matrices 
can  be  used  to  speed  implementation  by  coding  the  subroutines 
which  utilize  them  in  assembly  language.  The  difficulty  in 
deriving  the  scatter  terms  may  become  an  insurmountable  ob¬ 
stacle  when  non-isotropic  scatter  is  considered.  The  usual 
approach  used  for  non-isotropic  scatter  is  to  expand  the  flux 
in  Legendre  polynomials.  This  will  raise  the  order  of  the 
terms  in  the  scatter  integrals  and  significantly  increase  the 
difficulty  of  evaluating  them.  Since  the  local  matrices  are 
needed  for  the  penalty  calculation,  it  is  necessary  either  to 
save  them  or  to  recompute  them.  The  former  is  memory-inten¬ 
sive;  the  latter  increases  the  execution  time.  The  process 


of  making  a  mesh  by  hand  is  time-consuming  and  error-prone. 
The  lack  of  automatic  mesh  refinement  means  that  large  meshes 
are  not  readily  attainable.  Lastly,  the  columnar  scheme 
reduces  the  flexibility  of  the  finite  element  method,  and 
thus  increases  the  required  number  of  numerical  computations 
needed  for  a  given  level  of  accuracy.  It  does  not  seem  to  be 
a  severe  limitation,  but  specialized  applications  might  find 
it  cumbersome. 

There  are  many  benefits  to  be  gained  from  using  finite 
elements  to  solve  the  transport  equation:  the  ability  to 
concentrate  the  numerical  calculations  in  regions  where  the 
solution  is  rapidly  changing,  the  symmetry  and  sparseness  of 
the  global  matrix,  and  the  speed  of  the  calculations  involved 
in  the  finite  element  approximation.  This  study  has  charted 
a  path  toward  the  use  of  the  finite  element  method  that  hope¬ 
fully  will  be  of  benefit  to  those  who  may  choose  to  advance 
the  technique  even  further. 

Recommendations 

The  three  potential  sources  of  the  error  in  the  solution 
should  be  investigated  first.  Once  a  working  program  is 
achieved,  the  way  is  then  clear  to  improve  and  expand  upon 
the  technique  explored  in  this  study. 

The  basic  technique  explored  in  this  study  could  be 
improved  in  two  ways.  First,  the  amount  of  mathematical 


manipulations  needed  to  derive  the  finite  element  equations 
for  the  scatter  terms  should  be  reduced.  This  can  be  done  by 
giving  some  consideration  to  the  development  of  theorems  or 
formulas,  such  as  the  triangular  integration  formula.  For 
the  Legendre  polynomial  expansion  o  '  the  flux,  these  formulas 
may  indeed  by  a  necessity.  Additionally,  the  use  of  such 
theorems  should  increase  the  confidence  in  the  results  of  the 
derivation . 

Second,  alternatives  to  the  columnar  scheme  should  be 
considered.  This  does  not  seem  to  be  a  high  priority,  but 
the  current  scheme  does  somewhat  reduce  the  relative  advan¬ 
tage  of  the  finite  element  method  because  it  reduces  its 
ability  to  concentrate  numerical  calculations  in  selected 
regions. 

There  are  a  few  further  investigations  that  can  be  done 
to  increase  the  range  of  problems  which  can  be  handled.  Two 
immediate  needs  are  the  inclusion  of  the  finite  element  equa¬ 
tions  for  sources,  and  the  derivation  of  the  finite  element 
equations  for  non-isotropic  scatter.  This  method  should  be 
considered  incompletely  tested  until  an  automatic  mesh 
refinement  scheme  is  developed.  With  this  addition,  this 
method  could  be  tested  quantitatively  for  speed  and  accuracy. 

The  ultimate  goal  would  be  the  extension  to  two  dimen¬ 
sions.  Occasionally,  a  numerical  technique  works  well  in  one 
dimension,  but,  for  one  reason  or  another,  cannot  be 
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This  appendix  contains  the  derivations  of  the  five 


intervals  presented  in  Table  I  in  Chapter  III. 


The  first  integral  is 


Jli  <14  = 


where  /\  is  the  area  of  the  local  element.  Expressing 
the  integrals  indexed  by  l  as  a  vector,  yields 


or,  for  compatibility  with  later  integrals,  as  a  matrix 
tor  product,  and  with  different  scaling, 
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where  JT  is  the  unit  vector. 

The  second  integral  becomes 
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which,  in  vector  form,  is 
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which,  from  the  triangular  integration  formula  with  scaling, 
is  simply 
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which,  from  previous  definitions,  is 
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where 
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The  third  integral,  by  symmetry  with  the  second,  is 
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where 
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The  fourth  integral  becomes 
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from  the  triangular  integration  formula 
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In  matrix  form,  where  the  matrix  indices  are  over  j  and  )c 
not  over  L , 
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Including  all  three  matrices,  the  fourth  integral  is 


Explicitly  bringing  out  the  columnar  nature  of  the  integral 
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where  the  first  row  corresponds  to  (  =  1,  the  second  row  cor¬ 
responds  to  t  =  2,  etc.  Changing  the  form  of  the  expression 
by  multiplying  the  first  vector  with  the  matrix,  then  re¬ 
arranging  and  combining  terms,  yields 
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but  since  V2  =  ^3  »  then 
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fifth  term,  by  symmetry  with  the  fourth,  is 
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simply 
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PROGRAM  SATE 

* 

*  SOLUTION  OF  A  SELF-ADJOINT  TRANSPORT  EQUATION  BY  FINITE  ELEMENTS 

* 

IMPLICIT  UNDEFINED(A-Z) 

C  tt  Indexes. 

I NTEGER  I , J , K , TRI ANG , NONLOC , LOCMAT 
C  #*  Output  flags. 

LOGICAL  FMESH,FNODES,FELEM,FMTRXI ,FHTRXB»FFORCE 

C  **  Mesh  parameters. 

INTEGER  MINOD, MBNOD, MNTRIA, MLOCMT, SIZE 

C  t%  Maximum  dimensions  of  the  global  matrix  and  number  of  elements. 


PARAMETER  ( 

MINOD 

= 

64 

+ 

MBNOD 

S 

16 

+ 

MNTRIA 

64 

+ 

MLOCMT 

= 

512 

+ 

SIZE 

c 

MINOD+MBNOD 

C  **  Sise  data  st 

rue tu res. 

I NTEGER  PTNODE <  MNTRI A , 4 ) , COLUMN (32,2) 

REAL  C0RDND(MIN0D+MBN0D,2> ,PSI (MBNOD) ,PHI (MINOD) 

REAL  AREAS(MNTRIA) ,PENLTY (MNTRIA) 

REAL  M 1 1 ( MINOD , MI NOD ) , M i 2 ( M I NOD , MBNOD ) 

REAL  M21 ( MBNOD , MI NOD ) , M22 ( MBNOD , MBNOD ) 

REAL  MAT(MLOCMT ,3,3) 

C  ##  Mesh  variables. 

INTEGER  INTNOD,BNDNOD,NTRIAN 
REAL  RANGE, SIGMAT, SIGMAS 

C  tt  Computations. 

INTEGER  IER 
REAL  AREA,M(3,3> 

DOUBLE  PRECISION  MM(MINOD, MINOD) ,B(MINOD> ,UK(3*MIN0D> 

LOGICAL  BUG 

%%%%%%%%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%%%%*% 


C  #*  Initialise  the  element  characteristics. 

CALL  I NITAL ( MNTRI A , M I  NOD , MBNOD , NTRI AN , I NTNOD , BNDNOD , 
+  PTNODE, COLUMN, CORDND, PHI, PSI, AREAS, PENLTY, 

+  M11,M12,M21,M22, RANGE, SIGMAT, SIGMAS, B) 

C  **  Setup  the  output  flags. 

FMESH  =  .TRUE. 

FNODES  =  .TRUE. 


•*  *•  * 


FMTRXI  =  .TRUE 
FMTRXB  =  .TRUE 
FFORCE  =  .TRUE 


C  **  Assemble  the  matrix. 

LOCMAT  =  1 

DO  300  TRIANG=1,NTRIAN 
AREA  =  AREAS (TRIANG) 

DO  310  1=1,3 
DO  310  J=1 ,3 

MAT (LOCMAT , I , J)  =  0. 

310  CONTINUE 

CALL  STREAM  <  TR I ANG ,  AREA ,  MNTR I A ,  S I ZE ,  F'TNODE ,  CORDND ,  M ) 

CALL  ASEMBL ( TRIANG , TRI ANG, MNTRI A , MINOD , MBNOD ,  INTNOD , BNDNOD , 

+  PTN0DE,M,M11 ,M12,M21 ,M22) 

DO  320  1=1,3 
DO  320  J=l,3 

320  MAT (LOCMAT, I, J)  =  MAT (LOCMAT, I ,J)  +  M(I,J) 

CALL  ABSORB ( TRI ANG , ARE A , S ! GMAT , M ) 

CALL  ASEMBL ( TRIANG , TRIANG . MNTR I A , MI NOD , MBNOD , INTNOD , BNDNOD , 

+  PTN0DE,M,M11,M12,M21,M22> 

DO  330  1=1,3 
DO  330  J=l,3 

330  MAT (LOCMAT, I, J)  =  MAT (LOCMAT, I ,J)  +  M(I,J> 

LOCMAT  =  LOCMAT  +1 

NONLOC  =  COLUMN (PTNODE (TRI ANG, 4 ),1) 

DO  355  K=1,C0LUMN(PTN0DE(TRIANG,4) ,2) 

CALL  SCAT1 ( TRI ANG , NONLOC , MNTRI A , SI ZE , PTNODE , CORDND , ARE AS , 

+  SIGMAT, SIGMAS, M) 

CALL  ASEMBL ( TRIANG , NONLOC , MNTRI A , MI NOD , MBNOD , INTNOD , BNDNOD , 
+  PTNODE, M, Mil, M12,M21,M22) 

DO  340  1=1,3 
DO  340  J=1 ,3 

340  MAT (LOCMAT, I, J)  =  MAT(L0CMAT,I , J)  +  M(I,J) 

CALL  SCAT2(TRIANG, NONLOC, MNTRIA, SIZE, PTNODE, CORDND, AREAS, 

+  SIGMAS, M) 

CALL  ASEMBL ( TRIANG , NONLOC , MNTRIA , M I NOD , MBNOD , INTNOD , BNDNOD , 
+  PTNODE, M, Mil, M12,M21,M22) 

DO  350  1=1,3 
DO  350  J=1 ,3 

350  MAT(LOCMAT,I,J>  =  MAT( LOCMAT ,1, J)  +  M(I,J> 

CALL  SCAT3 ( TRIANG , NONLOC , MNTRIA , SIZE , PTNODE , CORDND , AREAS , 

*  SIGMAS, M) 
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CALL  ASEMBL ( TRIANG , NONLOC , MNTRI A , MINOD , MBNGD , INTNOD , BNDNOD 
+  PTNODE»M,Hll ,M12»M21 ,M22) 

DO  360  1*1,3 
DO  360  J=1 ,3 

360  MAT (LOCMAT , I , J )  =  MATCLOCMAT, I , J)  +  M(I,J> 


NONt.OC  =  NONLOC  +1 
LOCHAT  =  LOCMAT  +1 
355  CONTINUE 

C  —  Record  the  assembling  of  the  matrix. 

CALL  DEBUG ( 'SATE  * ,BUG> 

IF (BUG)  THEN 

PRINT  *,  ‘SATE***  TRIANG=  * » TRIANG 

CALL  OUTPUT <  MNTRI A , MINOD , MBNOD , NTRIAN , INTNOD , BNDNOD , PTNODE , 
+  CORDND, PHI ,PSI, AREAS, FENLTY, Mil, M12,M21 ,M22,B, 

+  .FALSE. , .FALSE. , .FALSE. ,FMTRXI , .FALSE. , .FALSE. ) 

ENDIF 

300  CONTINUE 


C  **  Compute  the  force  vector. 

DO  400  1=1, INTNOD 
P(I)  =  0. 

DO  400  J=l, BNDNOD 

B(I>  =  B<I)  -  M12<I,J)*PSI(J> 
400  CONTINUE 


C  **  Record  the  pre-solution  status. 

CALL  DEBUGCSATE  *,BUG> 

IF<BUG)  THEN 

CALL  OUTPUT(MNTRIA, MINOD, MBNOD, NTRIAN, INTNOD, BNDNOD, PTNODE, 
+  CORDND, PHI, PSI, AREAS, PENLTY, Ml 1 ,M12,M21 ,M22,B, 

+  FMESH ,F NODES ,FELEM , FMTRXI , FMTRXB , FFORCE ) 

ENDIF 


C  **  Convert  to  double  precision  for  IMSL.  Also,  advert  decompositi 
DO  500  1=1, INTNOD 
DO  500  J*l, INTNOD 
MM( I , J)  =  Mil (T,J) 

500  CONTINUE 


C  **  Solve  the  SATE  matrix  equation. 

CALL  LEQIF (MM, MINOD, INTNOD, INTNOD, B, MINOD, 1 ,0,WK, I ER) 
WRITER, '  (/,LEQIF,/'IER=  ',17,/)')  IER 

C  ##  Convert  the  LEQIF  returned  nodal  values  from  B  to  PHI. 
DO  600  1=1, INTNOD 
PHI < I )  =  B< I ) 

600  CONTINUE 


C  %%  Compute  the  penalty  of  each  element 


CALL  PEN  AL  ( NTR I  AN  >  MNTR I A ,  M I  NOD ,  MRNOD ,  MLOCMT ,  I  NT  NOIi , DNDNOD , PT NODE , 

+  COLUMN,  PHI,  PSI,  HAT,  Ml  1 ,M12,M21 ,M22,PENLTY> 

C  tt  Output  the  results. 

CALL  OUTPUT ( HNTRI A , MI NOD , MBNOD , NTRI AN , INTNOD , DNDNOD , PTNODE , 

+  CORDNU , PHI , PS I , AREAS , PENLTY , Ml 1 , M 12 , M2 1 , M22 , B , 

+  FMESH,FNQDES,FELEM,FMTRXI ,FMTRXF,FFORCE) 

jmm*mmmmmmmmmmmmrm**m***********m*mm 

STOP 

END 
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SUBROUTINE  DOCUM 

tt***X*******t%*%**t*******M**X****tttt*t*t*%*tt%*X%*X****t*t*t********* 
*  * 

*  DOCUMENTATION  * 

*  * 

ttt*%tt%t*1/;Wtt*tt*ttt%t*M*#**t****%M*tt*XVW*t**'X*%X*t*t******t*tt*lf-W 


lLt  Allan  Goff 
December  31,  1983 


* 

* 

* 

* 

He 

*  This  program  solves  the  self-ad.joint  transport  equation  in 

*  a  finite  element  scheme. 

* 


* 

* 

* 

* 

* 

* 

* 

* 


ttt**t****tt********t**$*********t**********x**************************** 

ttttttttttt************************************************************** 


*  SUBROUTINE  GLOSSARY  •  * 

*  * 

Ht  SATE  -  The  colling  program.  * 

*  DOCUM  -  The  documentation  routine.  He 

*  INITAL  -  Initializes  data  structures  and  element  constants.  * 

*  STREAM  -  Computes  the  contribution  fro®  the  streaming  term.  * 

*  ABSORB  -  Computes  the  contribution  from  the  absorption  term.  * 

*  SCAT1  -  Computes  the  contribution  from  the  1st  scatter  term,  * 

He  SCAT2  -  Computes  the  contribution  from  the  2nd  scatter  term.  * 

*  SCAT3  -  Computes  the  contribution  from  the  3rd  scatter  term.  * 

*  CASEDT  -  Determines  symmetry  or  assymmetry  between  elements.  He 

Ht  ASEMBL  -  Assembles  the  global  matrix  from  the  local  matrices.  * 

*  PENAL  -  Computes  the  minimization  integral  for  each  element.  # 

*  OUTPUT  -  Prints  selected  intermediate  and  final  data.  * 

*  DEBUG  -  A  low  over-head  debugging  routine.  * 

*  * 
*%tttttt*tt4t%**M****tttt***%***t%**M****t*****%****t******M%*ttt***** 

Htft****mmmmmm*mmnm*mmm**********m***********t* 


* 

* 

He 

* 

* 

* 

Ht 

Ht 


PARAMETER  GLOSSARY 

MNTRIA  -  Maximum  number  of  elements. 

MINOD  -  Maximum  number  of  interior  nodes. 

MBNOD  -  Maximum  number  of  boundary  nodes. 

MLOCMT  -  Maximum  number  of  local  matrices. 

PI  -  Ratio  of  diameter  of  a  circle  to  its  circumference. 


* 

* 

* 

Ht 

* 

Ht 

* 

* 


#$$t$$ttM#)iiW)nii[****#**tt$*t*$***************t*tt***t*tt********t******t 


$tt*t**tt$*t**t*t***t*tt*ttt*tttt*tt*$t*tt**t*ttttt$**t**t*t*t*tt*tt$ttt# 

*  VARIABLE  GLOSSARY  * 

*  Ht 

*  NTRIAN  -  Number  of  elements  in  this  mesh.  Ht 
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*  INTNOD  - 

*  BNDNOD  - 

*  LOCMAT  - 

*  LOC 

*  SJGMAT 

*  SIGMAS 

*  RANGE 

*  I ,  J,K  - 

*  II 

*  JJ 

*  TRIANG  - 

*  NONLOC  - 

*  X<3) 

*  U<3> 

*  A<3)  - 

*  B  ( 3  > 

*  C(3) 

*  AREA  - 

*  AREAS  O- 

*  PHIO  - 

*  PSIO  - 

*  PHIL0C(3) 

*  PHIPRM(3) 

*  Mil  ( , )  - 

*  M12< , )  - 

*  M21(,)  - 

*  M22( , 1  - 

*  MM<,>  - 

*  M(3,3)  - 
»  MAT  < , » )- 

*  BO 

*  Ml  ( ,  >  - 

*  MA(,>  - 

*  MB( , )  - 

*  MC( f )  - 

*  MD(,)  - 

*  PTN0DE(,3) 

*  CORDNIK  r2) 

*  COLUMN ( »2) 

*  SIZE  - 

*  P 

*  PP 

*  PEN 

*  PENLTYO 

*  UKO 

*  CASE  - 

*  BUG 

*  FMESH  - 

*  FNODES  - 

*  FELEM  - 

*  FMTRXI  - 

*  FMTRXB  - 


Number  of  interior  nodes  in  this  mesh.  # 
Number  of  boundary  nodes  in  this  mesh.  # 
Number  of  local  matrices.  * 
Number  of  local  matrices.  * 

-  The  total  macroscropic  cross  section.  # 

-  The  absorption  macroscropic  cross  section.  * 

-  Range  of  the  x-axis  in  units  of  sigmat.  # 
Indexes  for  general  use.  * 
The  row  index  on  the  global  matrix.  * 
The  column  index  on  the  global  matrix.  # 
Integer  specifying  the  current  local  triangle.  * 
Integer  specifying  the  current  non-local  triangle.  * 
X-coordinate  of  the  ith  vertex  of  the  current  element.  * 
Mu-coordinate  of  the  ith  vertex  of  the  current  element.  # 
X  coefficient  in  the  linear  finite  element  approximation.  * 
Mu  coefficient  in  the  linear  finite  element  approximation.  * 
Constant  coefficient  in  the  linear  finite  element  approx.  * 
Area  of  the  current  triangle.  * 
Areas  of  the  triangles.  t 
Vector  of  the  global  interior  nodal  values.  * 
Vector  of  the  global  boundary  nodal  values.  % 

-  Vector  of  the  local  interior  nodal  values.  * 

-  Vector  of  the  non-local  interior  nodal  values.  t 
The  global  matrix*  interior  node  section.  # 
The  global  matrix*  boundary  node  section,  off  diagonal.  * 
The  global  matrix*  boundary  node  section,  off  diagonal.  # 
The  global  matrix,  boundary  node  section,  on  diagonal.  * 
The  global  matrix*  interior  node  section,  double-p  copy.  * 
The  local  matrix.  * 
Record  of  local  matrices.  Used  to  compute  the  penalty.  * 
Force  vector  of  the  SATE  matrix  equation.  t 
Integer  matrix  from  scatter  terms,  X(1 )-coeff icient.  # 
Integer  matrix  from  scatter  terms,  X( 1 ) SX(2)-coeff icient.  * 
Integer  matrix  from  scatter  terms,  constant  coefficient.  # 
Integer  matrix  from  first  scatter  term,  symmetric  cases.  * 
Integer  matrix  from  first  scatter  terms,  antisymmetric  * .  * 


-  Pointers  to  the  global  nodes  from  the  element  vertexes.  * 


-  Coordintes  of  the  global  nodes. 

-  Number  of  elements  in  each  column. 

The  total  number  of  global  nodes. 

Common  foctors  in  the  local  matrix. 

Common  factors  in  the  local  matrix. 

The  penalty  value  from  the  worKing  element. 

-  The  penalty  values  for  all  the  elements. 

The  3*PIN0D  worKing  space  for  the  IMSL  routine  LEQIF. 
Indicates  symmetric  or  assymetric  case  between  elements. 
Debugging  flag  used  to  select  the  debugging  routine. 
Output  option  flag  for  mesh  definition. 

Output  option  flag  for  the  nodal  vectors. 

Output  option  flag  for  element  charactoristics. 

Output  option  flag  for  global  matrix. 

Output  option  flag  for  global  matrix. 


.  *k.  1 


*1 
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SUBROUTINE  INITAL(MNTRIA,MINCD,MBNOD,NTRIAN,INTNOD, BNDNOD, 

+  PTNOBE ,  COLUMN  , CORDND .  PHI ,  F3I ,  AREAS, PENLTY , 

+  M11,M12.M21,M22, RANGE, SIGMAT, SIGMAS, B> 

*  * 

*  INITIAL  COMPUTATIONAL  DATA 

* 

IMPLICIT  UNDEFINED(A-Z) 

C  **  Indexes, 

INTEGER  I,J,K,TRIANG,NODE 

« 

C  **  Passed  variables, 

INTEGER  MNTR I A , M I NOD , MBNOD 

INTEGER  NTRI AN , INTNOD , BNDNOD , NCOL 
INTEGER  PTN0DE(MNTRIA,4) ,CQLUMN(32,2) 

REAL  CORDND ( MI NQDiMBNOD , 2 ) 

REAL  PHKMINOD)  ,PSI  (MBNOIO 
REAL  AREAS(MNTRIA) , PENLTY (MNTRIA) 

REAL  Mil  (MIN0D,MIN0D),M12(MIN0Ii,  MBNOD) 

REAL  M21 (MBNOD, MINOD) ,M22( MBNOD, MBNOD) 

REAL  RANGE, SIGMAT, SIGMAS 
DOUBLE  PRECISION  B(MINOD) 


C  **  Computations, 

REAL  AREA,U(3) ,X(3) 

CHARACTER  TRASH116 
LOGICAL  BUG 


^  ^  ^  V  ^  ^  ^  X  ^  d-  ^  ^  V  dr  ^  ^  d  d  d  -u  d  d  d  d  d  d  d  d  d  d  ^  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  dr  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d  d 

^  ^  ^  T,  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  *  *  *  *  *  *  *  ^  *  *  *  *  *  ^  ^  ^  ^  ^  ^  ®  *•  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 


C  **  Read  in  the  mesh  definition, 

OPEN ( 3, F II E=* Data/mesh ' ) 

REWIND  3 

REAIK3, ' < A16) ' >  TRASH 

RE AD ( 3 , ' ( 4  < 1 X , 1 7  > ) ' )  NTR I AN , I NTNOD , BNDNOD , NCOL 
READ (3, ' (IX) ' ) 

READ(3, ' (A16) ' )  TRASH 

READ (3, ' (3< 1X,F7.3) ) ' )  RANGE , SI GMAT , SI GMAS 
RANGE  =  RANGERS I GMAT 
READ(3, '  (IX) ' ) 

READ(3, ' (A1A) ' )  TRASH 
DO  50  1=1 ,NTRIAN 

F;EAD(3, '(1X,I7,8X,4(1X,I7))')  TRIANG, (PTNODE( I , J) , J=1 ,4) 
50  CONTINUE 

READ(3, '  < IX) ' ) 

.V*.* 


92 


■*  * 


READ(3, ' < A16) ' )  TRASH 
DO  60  1=1 ,NCQL 

READ (3, ' <3( IX, I7,8X) ) ' )  TRIANG, (COLUMN( I , J) , J=1 ,2) 

60  CONTINUE 

READ(3,'(1X)'> 

READ(3, ' (A16) ' )  TRASH 
DO  70  1=1 , INTNOD+BNDNOD 

READ (3, ' (1X,I7,8X,2( 1X,F7.3) ) ' )  NODE, (CQRDND( I , J) , J=1 ,2) 
C0RDND(I,1)  =  CORDND (1,1) GRANGE 
70  CONTINUE 

READ(3,'(1X>'> 

REAIK3, ' < A16) 7  >  TRASH 
DO  80  1=1 , BNDNOD 

READ (3, ' (1X»I7,8X,1 (1X»F7.3) ) ' )  NODE,PSI(I> 

80  CONTINUE 

READ(3,'<1X)') 

CLOSE (3) 


C  **  Compute  the  ureas  of  each  of  the  triangular  elements. 

DO  100  TRIANG=1 ,NTRIAN 

U(3)  =  CORDND (PTNQDE( TRIANG, 3) ,2) 

U(2)  =  C0RDND(PTN0DE(TRIANG,2>,2) 

X (2)  =  CORDND ( PTNODE( TRIANG, 2 ),1) 

X(l)  =  CORDND(PTNODE< TRIANG, 1),1) 

AREA  =  ABS(  ,5*<  U<3)-U(2)  >*<  X(2)-X(l)  )  ) 

IF (AREA  ,LT.  1.0E-15)  AREA  =  1.0E-15 
AREAS(TRIANG)  =  AREA 
F’ENLTY(  TRIANG)  =  0. 

100  CONTINUE 

C  **  Zero  the  global  matrix  and  the  nodal  vector, 

DO  200  1=1 , INTNOD 
PHI ( I )  =  0. 

DO  21 0  J=l, INTNOD 
210  Mil ( I , J)  =  0. 

DO  220  K=1 ,BNDNOD 
M12( I ,K>  =  0. 

220  M21<K»I>  =  0. 

200  CONTINUE 

DO  230  1=1 ,BNDNOD 
DO  230  J=1 ,BNDNOD 
M22( I , J>  =  0. 

230  CONTINUE 

************************************************************************ 
CALL  DEBUG ( * INITAL* ,BUG) 

IF (BUG)  THEN 

PRINT  *,  * INITAL**** 

CALL  OUTPUT (MNTRI A, MI NOD, MRNQD, NTRI AN, INTNOD, BNDNOD, PTNODE, 


*  «  Jf 


SUBROUTINE  ASEMBL(TRIANG, NONLOC, MNTRI A, MINOD, MBNOD, INTNOD, BNDNOD, 

+  PTNODEf M»M11 ,M12,M21 ,M22) 

******************************************************************4****** 

* 

Assemble  the  Global  Matrix  * 

* 

j**************************** ****** ************ ****** ********** *********** 
IMPLICIT  UNDEFINED(A-Z) 

INTEGER  TR I ANG  , NONLOC , MNTRI A , MINOD  , MBNOD , INTNOD  , BNDNOB 
INTEGER  PTN0DE(MNTRIA,4) 

REAL  M( 3,3) 

REAL  Mil (M1N0B»MIN0D)  ,M12< MINOD, MBNOD) 

REAL  M21 ( MBNOD , MINOD ) , M22 ( MBNOD , MBNOD ) 

INTEGER  I,J,II,JJ 
LOGICAL  BUG 

************************************************************************* 


C  **  Relocate  the  local  matrix  elements  into  the  global  matrix* 

DO  100  1=1,3 
DO  100  J=1 , 3 

II  =  PTNODE<TRIANG, I ) 

JJ  =  PTNODE< NONLOC, J) 

IF<  (II  .LE.  INTNOD)  .  ANJ.  (JJ  .LE,  INTNOD)  )  THEN 
M1KII ,  JJ)  =  M11(II,JJ>  +  M(  I,  J) 

ELSE  IF(  (II  .LE.  INTNOD)  .AND,  (JJ  .GT.  INTNOD)  )  THEN 
M12< II , JJ-INTNOD)  =  M12( II , JJ-INTNOD)  T  M(I,J) 

ELSE  IF (  (II  .GT.  INTNOD)  .AND.  (JJ  .LE.  INTNOD)  )  THEN 
M21( II-INTNOD, JJ)  =  M21 ( I I-INTNOD, JJ)  +  M(I,J> 

ELSE  IF(  (II  .GT.  INTNOD)  .AND.  (JJ  .GT.  INTNOD)  )  THEN 
M22UI-INTN0D,  JJ-INTNOD)  =  M22(  II-INTNOD,  JJ-INTNOD) 

+  +M(I , J) 

ELSE 

PRINT  *,  'Impossible  option  in  ASEMBL.' 

ENDIF 

100  CONTINUE 

************************************************************************* 
CALL  DEBUG C ASEMBL', BUG) 

IF(BUG>  THEN 

PRINT  *,  'ASEMBL***' 

DO  8800  11=1, INTNOD+BNDN0D 
IF(II  .LE.  INTNOD)  THEN 

WRITE!*, 6210)  (Ml  1 < II , J) ,J=1, INTNOD ) , (M12( II , J) , J=1 ,BNDNOD) 
6210  FORMAT (1X,'!',20(1X,F6,3)) 

WRITE!*, 6220) 

6220  FORMAT! IX, ' ! ' > 

ELSE  < 

I  =  II  -  INTNOD 


SUBROUTINE  STREAM(TRIANG, AREA, MNTRIA, SIZE, PTNOD£.CORDND,M> 

*  * 

*  STREAMING  TERM  * 

*  * 
************************************************************************* 

IMPLICIT  UNDEFINED! A-Z) 

INTEGER  TRIANG,MNTRIA , SIZE ,PTN0DE(MNTRIA,4 ) 

REAL  AREA,C0RDND(SIZE,2) 

INTEGER  I,J 

REAL  U(3),A(3),P,PP,M<3,3) 

LOGICAL  BUG 

************************************************************************* 
C  **  Compute  the  mu-coordinates  of  eoch  vertex. 

U(l)  =  C0RDND(PTN0DE!TRIANG,1>,2) 

U(2)  =  CORDND ( FTNODE !  TRI ANG , 2 ) , 2 ) 

U(3)  =  C0RDND!PTN0DE!TRIANG,3),2) 

C  **  Compute  the  common  factors  in  the  local  matrix. 

P  =  U!1)*(U(1)+U<2>>  +  U(2)*(U(2)+U<3))  t  U!3)*(U(3)HI(  1) ) 

PP  =  SORT (P/( 24.* AREA)) 

C  **  Compute  the  x-coordinate  coefficient  for  each  vertex. 

A(l)  =  PP*(U(2)-U<3) ) 

A(2)  =  PP*(U!3)-U( 1 ) ) 

A(3)  =  PP*(U! 1 )-U(2) ) 


C  **  Compute  the  matrix  elements. 

DO  100  1=1,3 
DO  100  J=l,I 

M(I,J)  =  A(I)*A!J) 

M(  J,  I )  =  M<  I ,  J) 

100  CONTINUE 

************************************************************************* 
CALL  DEBUG ( *  STREAM 1 , BUG ) 

IF (BUG)  THEN 

PRINT  *,  'STREAM***  TRIANG=' , TRIANG 

DO  8800  1=1,3 

WRITE!*, 8600)  (M(I , J) , J=1 ,3) ,U( I ) , A(I) ,P,PP 
8600  FORMAT (1X,3(F7.3,1X),4X,4(F7.3,3X)) 

8800  CONTINUE 
ENDIF 
RETURN 
END 


SUBROUTINE  ABSORB<TRIANS,AREA,SIGMAT,M> 
************************************************************************* 

*  >K 

*  ABSORPTION  TERM  * 

%  * 

*********»*M****************************************  ******************* 

IMPLICIT  UNDEFINED(A-Z) 

INTEGER  TRIANG 
REAL  AREA, SIGMA T 

REAL  M(3,3> 

INTEGER  I,J 
REAL  PP 

LOGICAL  BUG 

###*#*#*##*###*######*#**##****###***###*####****#******#**#***#*##*#**## 


C  **  Compute  the  common  factors  in  the  local  matrix  elements* 
PP  =  SIGMAT*SIGMAT*AREA/12. 


C  **  Fill  in  the  local  matrix* 

Mil » 1 )  =  PP*2. 

M( 1 ,2)  =  PP*1. 

M<1,3>  =  PP*1. 

M<2,1>  *  PP*1. 
h<2,2)  =  PP*2, 

h(2,3)  =  PP *1. 

M(3, 1 )  =  PP*1. 

M(3r 2)  =  PP*1. 

M(3,3)  *  PP*2. 

************************************************************************* 
CALL  DEBUG  <’ ABSORB ', BUG ) 

IF (BUG)  THEN 

PRINT  <t,  'ABSORB###  TRIANG=* , TRIANG 
DO  8800  1=1,3 

WRITE**, 8600)  <M< I , J) , J=1 ,3) ,PP 
8600  FORMAT < 1X,3(F7 . 3,2X) ,5X, f 7,3) 

8800  CONTINUE 
ENDIF 
RETURN 
END 
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*  ■* 


SUBROUTINE  SCAT1 ( TRIANG ,NONLOC,MNTRIA, SIZE, PTNODE ,CORDND, AREAS , 

+  SIGMAT  .SIGMAS, M) 

t*tttt****xMt*tt**t**%tMtMt*****%*tt*t%*tt*t%%tt*t***%**t*tt*t***t*t** 

* 

First  Scatter  Terni  * 

*  * 
XX***X**tXttXXtXX**XtXXXtXX*Xt**t*X**t*XXX'i***X*tX*XXt*XXX***t*tX*XXtttt* 
IMPLICIT  UNDEFINED (A-Z) 


INTEGER  TRIANG, NONLOC,MNTRIA,PTNODE (MNTRIA,4) » SIZE 
REAL  AREAS(MNTRIA) ,CORDNIKSIZE, 2), SIGMAT, SIGMAS 

REAL  M<3,3> 

INTEGER  I, J, CASE 

REAL  MC<3,3>,MIK3,3),ALPHA,AREASQ,X<3>,D,P 
LOGICAL  BUG 

ttt******%***tt*******%**************tt******xt**x**t%***t*************** 
DATA  <MC(1,J),J=1,3>  /  4,  3,  3  / 

+  <MC<2,J),J=1,3)  /  3,  6,  6  / 

+  (MC(3, J) r J=lf 3)  /  3,  6,  6  / 

DATA  (MD(1,J),J=J,3)  /  4,  3,  3  / 

+  <MD(2,J),J=1,3>  /  3,  1,  1  / 

+  <HD<3,J),J>1,3)  /  3,  1,  1  / 

t*ttt**iii**tt*t***t*****tt**tt*t*t*t*ttt*t**tt*tt*tt**t%************%**tt* 

C  *#  Compute  the  coefficient  of  the  3rd  integral. 

ALPHA  =  SIGMASKSIGMAT  -  SIGMAS/2.) 


C  #*  Compute  the  product  of  the  areas  of  the  local  and  nonlocal  elements. 
AREASQ  =  AREAS ( TR I ANG ) *AREAS ( NONLOC ) 

C  ##  Compute  the  width  of  the  column  containing  these  elements. 

X(l)  =  C0RDND(PTN00£(TRIANG,1),!) 

X(2)  =  CORDNIKPTNODE (TRIANG, 2), 1) 

D  =  ABS(X<2)  -  X(l) > 


C  **  Compute  the  constant  factor  in  the  local  matrix. 
P  =  -ALPHA*AREASG/(30.*D> 


C  ##  Determine  the  case. 

CALL  CASEDT< TRIANG, NONL DC, MNTRIA, SIZE, PTNODE, CORDND, CASE) 

C  —  Symmetric  casef  fill  the  local  matrix  using  MB. 

IF1CASE  .EQ.  14)  THEN 
DO  100  1=1,3 
DO  100  J=l, 3 

H(I,J)  =  P*NC(I,J) 

100  CONTINUE 


C  —  Assymnietric  cose,  fill  the  local  matrix  using  Ml. 
ELSE  IF(CASE  .EG.  23)  THEN 
DO  110  1=1,3 
DO  110  J=1 ,3 

M(  I ,  J)  =  P*MD(I,J> 

110  CONTINUE 

ELSE 

PRINT  ' Inipossible  case  in  SCAT1.* 

ENDIF 


$tt*tM*tttttt********t*t****t*****%*tt*t**t*tt*****t*tt*t***t*****t***** 
CALL  DEBUG* *SCAT1  '  ,BU6) 

IF (BUG)  THEN 

PRINT  *,  *SCAT12»c**  triang=  * ,triang,*  NONLOC=',NONLOC 
PRINT  *,  'MtnJ)1 

WRITE**, ,(3(2X|F6.3))/)  <M( 1 , J) , J=1 ,3) 

WRITE**, '<3(2X,F6,3))')  <M(2, J) ,J=1,3) 

WRITE**, '(3(2X,F6.3))')  (M(3, J) , J=1 ,3) 

PRINT  *,  *D  P  ALPHA  AREASQ  * 

PRINT  *,  H,P, ALPHA, AREASQ 
ENDIF 


RETURN 

END 


SUBROUTINE  SCAT2( TRI ANG, NONLOC, MNTRI A, SIZE, PTNODE .CORDND, ARE AS, 

+  SIGMAS, M) 

*  * 

*  Second  Scatter  Tern  # 

#  * 

IMPLICIT  UNDEFINED(A-Z) 


INTEGER  TRI ANG , NONLOC , MNTR I A , PTNODE ( MNTRI A , 4 ) , SIZE 
REAL  AREAS (MNTRI A ),CORUNH< SIZE, 2), SIGMAS 

REAL  H<3,3> 

INTEGER  I, J, CASE 

REAL  MB <3, 3), Ml (3, 3) ,AJP(3) ,X(3),U(3),UP(3),B(3),D,P 
LOGICAL  BUG 

*tt*******t***M**tt****t*%tt**t*****ttt'll.%*t*tt*%t%*t*%%**t******tt*****t 
DATA  (MB( 1 , J) , J=1 ,3)  /  4,  3,  3  / 

+  <MB(2» J) , J=1 ,3)  /  3,  8,  4  / 

*  (MB(3, J) , J=1 ,3)  /  3,  4,  8  / 

DATA  (Ml (1 , J) , J~l,3)  /  6,  2,  2  / 

+  (Ml (2, J) , J=1 ,3)  /  2,  2,  1  / 

+  (Ml (3, J) , J=1 ,3)  /  2,  1,  2  / 

tt**ttttt*tMtt*tt%tt*tt*%t***%tt****tM*M***%t%tt**t*t*%*%*t**t*%***t*X 


C  **  Compute  the  width  of  the  column  containing  these  elements. 
X(l)  =  C0RDND( PTNODE (TRI ANG, 1),1) 

X(2)  =  C0RDND(PTN0BE(TRIANG,2> ,1) 

D  =  ABS(X(2)  -  X(l) ) 


C  **  Compute  the  mu  coordinates  of  the  local  element. 
U(l>  =  CORDND (PTNODE (TRI ANG, 1 >, 2) 

U(2)  =  CORDND (PTNODE (TRI ANG, 2), 2) 

U(3)  =  CORDND (PTNODE (TRI ANG, 3) ,2) 


C  *#  Compute  the  mu-prime  coordinates  of  the  non-local  element. 
UP(1)  =  CORDND ( PTNODE ( NONLOC, 1), 2) 

UP(2)  =  CORDND (PTNODE (NONLOC, 2), 2) 


UP(3)  =  CORDND (PTNODE (NONLOC, 3), 2) 


C  **  Compute  the  x  coefficient  of  the  non-local  element's  linear 
C  —  interpolating  function. 

AJP(l)  -  UP(2)  -  UP(3) 

AJP(2)  =  UP(3)  -  UP(1) 

AJP(3)  =  UP(1)  -  UP(2) 


C  tt  Compute  the  constant  factor  in  the  local  matrix. 
P  =  (SIGMAS/2.)  *  <AREAS(TRIANG)/(60,*D>) 


C  Determine  the  cose. 

CALL  CASEDT (TR1ANG, NONLOC, MNTRIA, SIZE, PTNODE,CORDND, CASE) 


C  —  Symmetric  case,  fill  the  locol  matrix  using  MB. 
IF(CASE  .EG.  14)  THEN 
DO  100  1=1,3 
DO  100  J=1 ,3 

100  B(I)  =  P*MB(I,J)*U(J) 

DO  110  1=1,3 
DO  110  J=1 *3 

110  M(I,J)  =  B<I)*AJP(J) 


C  —  Assymmetric  case,  fill  the  local  matrix  using  Ml. 
ELSE  IF (CASE  .EQ.  23)  THEN 
DO  120  1=1,3 
DO  120  J=1 ,3 

120  B(I)  =  P*M1<I,J)*U<J> 

DO  130  1=1,3 
DO  130  J=1 , 3 

130  M(I,J)  =  B< I )*AJP< J) 


ELSE 

PRINT  K,  'Impossible  case  in  SCAT2.* 

ENDIF 

**********m***m****»*m#*mm**#mmmm***mmm****m*** 

CALL  DEBUG CSCAT2  ’,8UG) 

IF(BUG)  THEN 

PRINT  *,  'SCAT2***  TRIANG=  ’ ,TRIANG, '  NONLOC=' , NONLOC 
PRINT  *,  '  M(I, J)  U  UP  AJP* 

WRITE (*, ' (6(2X,F£.3) ) '  >  (M(l , J) , J=1 ,3) ,U( 1 ) ,UP( 1 ) , AJP(1 ) 
WRITE(*,'(6(2X,F&.3))'>  <M(2, J) , J=I ,3) ,U(2) ,UP(2) ,AJP(2> 

WRITEU,  '<A(2X,F6.3>>'>  <M(3,  J) ,  J=1 ,3)  ,U(3)  ,UP(3> ,  AJP<3> 

PRINT  *,  ’D  P* 

PRINT  *,  D.P 
ENDIF 


RETURN 

END 


SUBROUTINE  SCAT3< TRIANG, NONLOC, MNTRIA, SIZE, PTNODE, CORDND, AREAS, 

+  SIGMAS, H) 

*  * 

%  Third  Scatter  Terra  * 

*  * 

IMPLICIT  UNDEFINED(A-Z) 

INTEGER  TRIANG,NONLOC,MNTRIA,PTNODE(MNTRIA, 4) ,SIZE 
REAL  AREAS<MNTRIA) ,CORBNIKSIZE, 2), SIGMAS 

REAL  M<3,3> 

INTEGER  I, J, CASE 

REAL  MC(3,3),MD(3»3)»X(3),UP(3), A JP (3>,B(3),D»P 
LOGICAL  BUG 

*%%*%%%%%%%%%%%*%%%%%%%%*%%%%%%*%*%%%%%%*%%%%%%%*%%%%%%%%%%%%%%%%%*%%%%%* 
DATA  (MC( 1 , J) , J=1 ,3)  /  4,  3,  3  / 

+  (MC(2, J) , J=1 ,3)  /  3,  6,  6  / 

f  (MC(3, J) , J=1 ,3)  /  3,  6,  6  / 

DATA  <MD(1 , J) , J=1 ,3)  /  4,  3,  3  / 
f  (MD(2» J) ,J=1,3)  /  3,  1,  1  / 

+  (MIK3,  J) ,  J=1 ,3)  /  3,  1,  1  / 

tt%ttt**tt***t%%*t*%*%*%**********t**t*****tt*t**t**tt*t**xt*******t**tt* 

C  ##  Compute  the  width  of  the  column  containing  these  elements. 

X(l)  =  CORDND ( PTNODE (TRIANG, 1 ) ,1 ) 

X<2>  =  CORDND (PTNODE (TRIANG, 2 >,1> 

D  =  ABS(X(2>  -  X(l)> 

C  tt  Compute  the  mu-prime  coordinates  of  the  non-local  element. 

UP(l)  =  CORDND ( PTNODE ( NONLOC , 1 ) , 2 ) 

UP ( 2 )  =  CORDND (PTNODE (NONLOC, 2), 2) 

UP(3)  =  CORDND (PTNODE (NONLOC, 3), 2) 

C  t*  Compute  the  >:  coefficient  of  the  non-local  element's  linear 
C  —  interpolating  function. 

AJP(l)  =  UP < 2 )  -  UP (3) 

AJP<2>  =  UP  ( 3 )  -  UP(1) 

AJP(3)  =  UP(1)  -  UP<2> 

C  tt  Compute  the  constant  factor  in  the  local  matrix. 

P  =  -(SIGMAS/2.)  *  <AREAS(TRIANG)/(ivO.*D)) 

C  #*  Determine  the  case. 

CALL  CASEDT < TRI ANG , NONLOC ,MNTRIA , SIZE, PTNODE , CORDND, CASE ) 


C  —  Symmetric  case,  fill  the  local  matrix  using  MB 
IF(CASE  .EQ.  34)  THEN 


no  100  1=1,3 
no  100  j=i ,3 

100  B<I)  =  P*MC(I,J>*UP(J> 

DO  110  1=1,3 
no  no  j=i ,3 

110  H(I,J>  =  B<I)*AJPCJ) 

C  —  Assynimetric  cose,  fill  the  local  matrix  using  Ml. 
ELSE  IF (CASE  .EG.  23)  THEN 
DO  120  1=1,3 
DO  120  J=1 ,3 

120  B<  I )  =  P*MD<I,J)#UP<J> 

DO  130  1=1,3 
DO  130  J=1 ,3 

130  H(I,J)  =  B(I)*AJP<J) 

ELSE 

PRINT  *,  .'Impossible  cose  in  SCAT3.' 

ENDIF 


:m**#*******#!m**!m*)m#:m!m:mm*:m:m*!m****#*#**#********#**#* 

CALL  DEBUG CSCAT3  ’.BUG) 

IF(BUG)  THEN 

PRINT  *,  'SCAT3***  TRIANG=  '  ,TRIANG, '  NONLOC=*  ,NONLOC 
PRINT  *,  *M( I , J)  UP  AJP' 

WRITE (*, '(5<2X»F6.3))')  <M< 1 , J) , J=1 ,3) ,UP< 1 ) , AJP( 1 ) 
URITE(*,/(5(2X,F6.3))/)  <M<2, J) , J=1 ,3) ,UP<2> ,AJPC2) 

WRITE  (<t,  ,(5<2X,F6.3))/)  ( M  ( 3 ,  J ) ,  J=  1 , 3 ) ,  UP  ( 3 ) ,  A  JP  ( 3 ) 

PRINT  *,  *D  P' 

PRINT  *,  D,P 
ENDIF 


RETURN 


SUBROUTINE:  CASEBT(TRIANG, NONLOC, MNTRIA, SIZE. PTNODE, CORDND, CASE) 

*#*  *****  *******  ******  ***********************************  ******** 
*  # 

*  Deterniine  case  of  the  local  and  non-local  elements.  * 

*  * 
************************************************************************* 

IMPLICIT  UNDEFINED (A-Z) 

INTEGER  TRI ANG ,  N0NL0C,MNTF;I  A ,  SIZE ,  F'TNQDE  ( MNTRIA ,  4 ) 

REAL  CORDND (SIZE, 2) 

INTEGER  CASE 

REAL  X ( 3 ) , Y ( 3  > , LOCL , NOLOCL 

LOGICAL  BUG 

************************************************************************* 


C  **  Compute  the  coordinates  of  the  columnar  ends  of  the  triangles. 

X(l>  =  CORDND(PTNOIE(TRIANG, 1 ) ,  1 ) 

X(2)  =  C0RDND(PTN0D£<TRIANG»2) , 1 ) 

Y(l)  =  CORDND ( PTNODE ( NONLOC ,  1 ) ,  1 ) 

Y<2)  =  CORDND (PTNODE (NONLOC, 2 ),1) 

C  **  Determine  the  direction  in  which  the  triangles  are  pointing. 

LOCL  =  X(2)  -  Xd) 

NOLOCL  *  Y(2)  -  Yd) 

C  **  Determine  the  case. 

IF(LOCL*NOLOCL  .GT.  0.)  THEN 
CASE  =14 
ELSE 

CASE  =  23 
ENDIF 

************************************************************************* 
CALL  DEBUG CCASEDT*, DUG) 

IF(BUG)  THEN 

PRINT  *,  ’CASEDT***  CASE=  ’,CASE 
PRINT  *,  *L0CL=  * ,L0CL, *  N0L0CL=  *, NOLOCL 

ENDIF 


RETURN 

END 
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SUBROUTINE  PENAL(NTRIAN,HNTRIA,MINGD,MBNOD.MLOCMT , INTNOD.DNDNOD, 

+  PTNODE, COLUMN, PHI, PSI, MAT, Ml 1,M12,M21,M22,PENLTY) 

#*#***  **##*##***#**#***#**************#  *'.***1**  *********  ****************** 

*  PENALTY  FUNCTION  * 

*  * 
************************************************************************* 

IMPLICIT  UNDEFINED(A-Z) 

INTEGER  NTRI AN , MNTRI A , MINQD , MBNOD , MLOCMT , I NTNOD , BNDNOD 
INTEGER  FTN0DE(MNTRIA,4>,C0LUMN<32,2> 

REAL  PHI(MINOD) , PSI (MBNOD) 

REAL  MAT (MLOCMT, 3, 3) 

REAL  Mil  <  MI  NOD, MI  NOD)  ,  M12(MIN0D, MBNOD) 

REAL  M21 ( MBNOD , MINOD  > , M22 ( MBNOD , MBNOD ) 

REAL  PENLTY(MNTRIA) 

INTEGER  I , J , K , TR I ANG , NONLOC , I I , J J , LOC 
REAL  PHI  LOC  (3)  ,  PHIF’RM(3)  ,F’EN 
LOGICAL  BUG 

I************************************************************************* 

CALL  DEBUG C PENAL  *,BUG) 

LOC  =  1 

C  ##  Compute  the  cummulotive  penalty  function  element  by  element* 

DO  100  TF;IAN6=1  »NTRIAN 
PEN  =  0. 

C  ##  Compute  the  local  terms  penalty 
C  —  Compute  this  element's  phi  vector. 

DO  110  1=1,3 

II  =  PTNODE<TRIANG, I ) 

IF < 1 1  *LE.  I NTNOD)  THEN 
PHILOC(I)  =  PHI (II) 

ELSE 

PHILOC(I)  =  PSI ( II-INTNOD) 

ENDIF 
110  CONTINUE 

IF (BUG)  THEN 

PRINT  *,  'PENAL#*#  TF;IANG=  MRIANG 

PRINT  *,  *  PEN  FHILOC(I)  MAT(LOC,I,J)  PHILOC(J)* 

ENDIF 

DO  120  1=1,3 
DO  120  J=1 ,3 

PEN  =  PEN  +  PHI LOC ( I )*MAT (LOC, I , J)#PHILOC<  J) 

IF (BUG)  THEN 

URITE(#,  6120)  F'EN,F'HILOC( I ) ,  I  ,MAT (LOC,  I ,  J) ,  J,F'HILOC(  J) 
FORMAT ( IX ,F7 *3, 7X , 3( F7. 3,2X , 13, 2X) ) 


6120 


120  CONTINUE 


C  ##  Compute  the  non-local  terms  penalty 
LOC  =  LOC  +  1 

NONLOC  =  COLUMN ( PTNODE ( TRIANG, 4 ) ,  1 ) 

DO  355  K=1 , COLUMN (PTNODE (TRIANG, 4)  ,2) 

C  —  Compute  this  element's  non-local  phi  vector. 

DO  130  J=l,3 

JJ  =  PTN0DE( NONLOC, J> 

IF( JJ  .LE.  INTNOD)  THEN 
PHIPRM(J)  =  PHI ( JJ) 

ELSE 

PHIPRM(J)  =  PSI ( JJ-INTNOD) 

ENDIF 

130  CONTINUE 

IF (BUG)  THEN 

PRINT  *,  "PENAL*##  TRIANG=  * ,TRIANG 

PRINT  #,  "  PEN  PHILOC(I)  MAT(L0C,I,J)  PHIPRM(J)" 

ENDIF 

DO  140  1=1,3 
DO  140  J-1,3 

PEN  =  PEN  +  PHIL.OC(I)#MAT(LOC,I,J)*PHIPRM(J) 

IF (BUG)  THEN 

WRITE(*,6120)  PEN,PHILOC(I),I,MAT(LOC,I,J),J,PHIPRr (J) 
ENDIF 

140  CONTINUE 

LOC  =  LOC  +  1 
NONLOC  =  NONLOC  +  1 
355  CONTINUE 

PENLTY(TRIANG)  =  ,5*PEN 
100  CONTINUE 


C  **  Compute  the  global  penalty. 

PEN  =  0. 

DO  200  11=1 , INTNQD+BNONOD 
DO  200  JJ=1 , INTNOD+BNDNOD 

IF (  (II  .LE.  INTNOD)  .AND.  (JJ  .LE.  INTNOD)  )  THEN 
PEN  =  PEN  +  PHI  ( II  )*H11  ( II » JJ)#F'HI(  JJ) 

ELSE  IF (  (II  .LE.  INTNOD)  .AND,  (JJ  ,GT .  INTNOD)  )  THEN 
PEN  =  PEN  +  PHI < II >*M12< II,  JJ-1NTN0P)*PSI ( JJ-INTNOD) 

ELSE  IF (  (II  .GT.  INTNOD)  .AND.  (JJ  .LE.  INTNOD)  )  THEN 
PEN  =  PEN  +  PSI  ( II-INTN0D)#M21  ( I I-INTNOD,  JJ)*F‘HI  ( JJ) 

ELSE  IF (  (II  .GT.  INTNOD)  .AND.  (JJ  .GT.  INTNOD)  )  THEN 
PEN  =  PEN  +  PSI ( I I-INTNOD )#M22( JJ-INTNOD, I I-INTNOD)# 

+  PSI (JJ-INTNOD) 


PRINT  *»  'Impossible  option  in  PENAL!  global  penalty.' 

ENDIF 

200  CONTINUE 

PEN  =  ,5*PEN 

*ttt.t*****1f-**t*^*^**%**********^*4.*****t*******%*****%***%**^*%il.*t*t*** 
IF (BUG)  THEN 

PRINT  *,  'PENAL  global  penalty  =  ',  PEN 
ENDIF 
RETURN 


SUBROUTINE  OUTPUT (MNTRIA ,MINOD. MBNOD, NTF;IAN,INTNOD,BNDNOD,PTNODE, 

+  CORDND, PHI, FSI, AREAS, F'ENLTY, Ml 1 ,M12,M21 ,M22,B, 

+  FMESH , FNODES , FELEM , FMTRXI , FMTRXB , FFORCE ) 

££££££££££££££££££££*££££££££££££*£££****££**£££££*£££*££££££££££££££££££ 
£  £ 

*  PRINT  SELECTED  ITEMS  * 

*  * 
££££££££££££X£££££££££££££££££££££££*£#££££***£*££££)|c#££££££*££££££££££££ 

IMPLICIT  UNDEFINED* A-Z) 

INTEGER  MNTRIA, M1NOD, MBNOD, NTRIAN,INTNOD,BNDNOD,PTNODE( MNTRIA, 4) 
REAL  AREAS(MNTRIA) , PENLTY* MNTRIA) ,CORDND(hINODtMBNOD, 2) 

REAL  PHI *MINOD) ,PSI (MBNOD) 

REAL  Mil (MINQD,MINOD) ,M12(MIN0D, MBNOD) 

DOUBLE  PRECISION  B*MINOD) 

REAL  M2 1 ( MBNOD , MINOD > , M22 ( MBNOD , MBNOD ) 

LOGICAL  FMESH , FNODES , FELEM, FMTRXI , FMTRXB , FFORCE 

INTEGER  I , J, K, TRIANG, I I 
REAL  PEN 
LOGICAL  BUG 

£££££££££££££££££££££££££*££££££££££££££*<£££*££*££££££££££££££££££££££££ 


£  There  are  six  groups  of  information.  Selected  by  flags.  * 
$  £ 
£  FNODES  -  Nodal  values.  £ 
£  FELEM  -  Element  characteristics.  £ 
£  FMTRXI  -  Elements  of  global  matrix,  (interior  part).  £ 
£  FMTRXB  -  Elements  of  global  matrix,  (boundary  part),  £ 
£  FFORCE  -  Force  vector  of  SATE  matrix  equation.  £ 
£  FMESH  -  Mesh  definition.  £ 
£  £ 


£££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££ 

C  ££  Mesh  definition. 

IF(FMESH)  THEN 

PRINT  £,  ’Mesh  Definition}  • 

PRINT  £,  *  Triangle  Global  Nodes.’ 

DO  50  TRIANG=1 ,NTRIAN 

WRITE(£,6050)  TRIANG, (PTNODE* TRIANG, I ), 1=1 ,4) 

6050  FORMAT (4X, I3,SX,4( I3,3X) ) 

50  CONTINUE 
PRINT  £ 

PRINT  £ ,  ’  Node  Coordinates  (x,mu)«* 

DO  60  1=1 , INTNODfBNDNOD 

WRIT E*£, 6060)  I , (CORDND* I , J) , J=1 ,2) 

6060  FORMAT ( 4X , 13 , 7X , 2 ( F7 » 3 » 3X ) > 

60  CONTINUE 
PRINT  £ 

ENDIF 


IF(FNODES)  THEN 

PRINT  *,  ‘Nodal  values.* 

DO  100  1=1 , INTNOD 

WRITE**, 6010)  I  ,F'HI (I ) 

6010  FORMAT (2X, 13 , 3X,F9.4 ) 

100  CONTINUE 

DO  110  1=1 , BNDNOD 

WRITE**, 6010)  I+INTNOIW F'SI  ( I ) 
110  CONTINUE 
PRINT  * 

ENDIF 


C  **  Element  characteristics. 

IF(FELEM)  THEN 

WRITE**, 6100) 

6100  FORMAT ( ‘Triangle  Area  Penalty  Nodes(x,u)‘) 

DO  200  TRIANG=1,NTRIAN 

WRITE**, 6110)  TRIANG, AREAS (TRIANG) f PENLTY(TRIANG) , 

+  ( (CORDND(PTNODE(TRIANG, J) ,K) ,K=1 ,2) , J=1 ,3) 

6110  FORMAT *1X,I3»3X,F6.3,3X,F7.4,4X,3*2(F6.3»1X),3X)) 

PEN  =  PEN  +  PENLTY( TRIANG) 

200  CONTINUE 

WRITE**, 6120)  PEN 
6120  FORMAT * 15X,F9.5) 

ENDIF 

C  **  Elements  of  global  matrix,  (interior  part). 

IF(FMTRXI)  THEN 
PRINT  it 
WRITE**, 6200) 

6200  FORMAT*"  Global  matrix,  (interior  part).*) 

PRINT  * 

DO  300  1=1, INTNOD 

WRITE**, 6210)  (M11(I,J),J=1, INTNOD) 

6210  FORMAT  <1X,‘!‘,20*1X,F6,3)) 

WRITE**, 6220) 

6220  FORMAT (IX, ' ! ‘ ) 

300  CONTINUE 
PRINT  * 

ENDIF 

C  **  Elements  of  global  matrix,  (boundary  part). 

IF(FMTRXEt)  THEN 

PRINT  it 

PRINT  it,  '  Global  matrix,  (boundary  part*  M12)." 

PRINT  * 

DO  400  1=1, INTNOD 

WRITE**, 6210)  ( M 1 2 ( I ,  J ) ,  J=  l , BNDNOD ) 

WRITE  (it,  6220) 

400  CONTINUE 


PRINT  * 

PRINT  *,  *  Global  matrix,  (boundary  part:  H21).' 
PRINT  * 

DO  410  1=1 ,BNDN0D 

WRITE ( *,6210)  (M21(I,J)»J=1, INTNOD) 

WRITE**, 6220) 

410  CONTINUE 

PRINT  * 

PRINT  *,  ’  Global  matrix,  (boundary  part:  M22).’ 
PRINT  * 

DO  420  1=1 ,BNDN0D 

WRITE* *,6210)  (M22(I,J),J=1,BNDN0D) 

WRITE**, 6220) 

420  CONTINUE 
PRINT  * 

ENDIF 


C  **  Force  vector  of  the  SATE  matrix  equation. 
IF(FFORCE)  THEN 

WRITER,' (’Force  vector.*)') 

DO  500  1=1, INTNOD 


WRITE<*,6500)  B( I ) 
6500  F0RMAT(2X,F7.3) 

500  CONTINUE 
PRINT  * 

ENDIF 


*  ^  ^  ^  V  v  T  ^  x  X  X  ^  X  X  X  X  ^  ^  X  X  X  ^  x  X  X  X  X  X  X  X  X  X  X  ^  ^  ^  ^  ^  ^  X  x  X  X  X  ^  ^  ^  X  X  ^  ^  ^  X  ^  X  ^  X  X  X  X  X  X  ^  ^  ^  ^  x  X  T  X  x 

CALL  DEBUG (’OUTPUT*, BUG) 

IF (BUG)  THEN 

DO  8800  11=1, INTNOD+BNDNOD 
IF(II  .LE.  INTNOD)  THEN 

WRITE<*,6210)  (Mil ( II, J),J=1, INTNOD), (H12( II, J) , J=1 ,BNDN0D) 
WRITE**, 6220) 

ELSE 

I  =  II  -  INTNOD 

WRITE**, 6210)  <M21(I,J),J=1, INTNOD),  ( M22 ( I , J ) , J= 1 , BNDNOD ) 
WRITE**, 6220) 

ENDIF 

8800  CONTINUE 
ENDIF 


RETURN 

END 


SUBROUTINE  DEBUG ( SUB NAM , BUG ) 

C  ****************************#*************************#**:****#*****#*#* 


C  ***  *#* 

C  ***  DOCUMENTATION  -  DEBUG  *** 

C  ***  *#* 

c  **#********************##*#****#*****#*##**#**#**#**#*#*****###******** 
c  *  * 

C  *  Lt »  Allan  Goff  August  27,  1983  * 

C  * - * 

C  *---  DEBUGGING  AID  * 

C  *— - * 

C  *  **  ABSTRACT  **  * 

C  *  * 

C  t  This  routine  scans  a  list  of  subroutine  names  to  be  debugged  * 

C  *  and  if  the  calling  routine  is  on  the  list  it  returns  a  true  value  * 

C  *  via  the  variable  BUG.  * 

C  * - * 

C  *  **  BACKGROUND  **  * 

C  *  * 

C  * - * 

C  *  ##  FLOW-CHART  **  « 

C  *  * 

C  * - - - * 

C  *  *«  IHF'LIMENTATION  **  * 

C  *  * 

C  *  The  list  of  subroutines  to  be  debugged  may  be  as  large  as  * 

C  *  100  names.  However,  DEBUG  quits  looKing  thru  the  list  when  it  * 

C  *  finds  the  name  QUITS8.  # 

C  * - * 

C  *  **  INTERFACE  REQUIREMENTS  **  * 

C  *  * 

C  t  The  parameters  passed  areJ  # 

C  *  SUBNAM  -  The  name  of  the  calling  program.  # 

C  *  BUG  A  The  flag  which  indicates  if  debuggins  is  to  be  done.  * 

C  * - y 

C  *  **  LOCAL  VARIABLES  **  * 

C  *  * 

C  t  BUGEM  -  The  list  of  subroutines  to  be  debugged.  # 

C  *  I  -  Index  variable  used  in  scanning  the  list.  % 

C  * - * 

C  *  *#  LIMITATIONS  **  *  * 

C  *  * 

C  * - » 

C  *  **  REFERENCES  **  * 

C  *  * 

C  ##*#*#****r**#***#***#*##**##***##**»*******»**#*****###*###******#**<* 

c  ***  *** 

C  ***  CODE  -  DEBUG  **# 

C  ***  *** 


c  *#*»mmmmmmmmmmmm#***#***#******#*****mmm* 

c  ##***m**mn**#*m*******  declarations  mm****mm******m*m 


C  -  HAVE  COMPILER  CHECK  FOR  MISSPELLED  VARIABLE  NAMES 

IMPLICIT  UNDEFINED(A-Z) 

c  -  INDEXES  - 

INTEGER  I 

C  -  PASSED  VARIABLES  - 

CHARACTER  SUBNAM*6 
LOGICAL  BUG 

C  -  LOCAL  VARIABLES  - 

CHARACTER  BUGEM<101)*6 
DATA  BUGEH<1)  /'GUIT88'/ 

DATA  BUGEM(2)  /'SATE  '/ 

DATA  BUGEM(3)  /'INITAL'/ 

DATA  BUGEMi 4)  /'ASEMBL'/ 


DATA  BUGEM(5)  /'STREAM'/ 

DATA  BUGEM(6)  /'ABSORB'/ 

DATA  BUGEH(7)  /'SCAT1  '/ 

DATA  BUGEM(8)  /'SCAT2  V 
DATA  BUGEMC9)  /'SCAT3  '/ 

DATA  BUGEM(IO)  /'CASEDT'/ 

DATA  BUGEM(ll)  /'PENAL  '/ 

DATA  BUGENU2)  /'OUTPUT'/ 

DATA  BUGEM( 13)  /'0UIT83'/ 

C  mmnimmm************  EXECUTABLES  t*tt*%X*ttttt*n************* 
BUG  =  .FALSE. 

1  =  0 

2000  I  =  I  F  1 

IF(FUG£M(I)  .EO.  '0UIT33')  GOTO  999 
IF •' BUGEM(  I )  .EG.  SUBNAM  )  THEN 
BUG  =  .TRUE, 

GOTO  999 
ENDIF 
GOTO  2000 

c  **m****#m*:m***:m*#*****#*  DONE  ;mm:mm:m*#*****m***#**#* 
999  RETURN 
END 
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Appendix  G 


In  this  appendix,  the  local  matrix  for  the  first  scat¬ 
ter  term,  for  the  antisymmetric  cases  (Cases  2  and  3),  is 
derived. 

The  starting  point  for  this  derivation  is  the  extremi- 
zation  integral. 

J3  =  ctu'dA  ^  (30i) 

Starting  with  the  usual  finite  element  approximation,  the 
flux  at  the  local  element  is 

*  5  * 

=  €  <Pi  =  £  f?  “  (*i  *  +  k  M  +  b  J  ( 302) 


regrouping  terms, 


(303) 


and  defining 


(304) 


(305) 


(306) 


c{tX  f  4,  A  f 


(307) 


Similarly,  the  finite  element  approximation  for  the  flux  at 
the  non-local  element  is 


—  A*  +A-*  +  Aa 


(308) 


where  the  coefficients  are  defined  as 

J 


(309) 


^  ^  L 


(310) 


*  =  § 


'  J 


(311) 


Thus,  the  first  scattering  term  extremization  integral 


becomes 


(312) 


where  and  are,  respectively,  the  upper  and  lower  edges 
of  the  local  element,  M*  and  44*  are,  respectively,  the  upper 
and  lower  edges  of  the  non-local  element,  and  . 
Integrating  over  M  and  M* 


I,  -  «/  f|[^(  +  X  rtl'H'J 


s 


[fo**  h)*'+ 


(313) 


147 


(314) 


(327) 


(328 


KX+ 

+  (4*»  *•-*  ~*i  «•**  f-)*,  +  «f3} 


+  i  ^/)  +■  *,*,  f 


(329 


let 


£,  =  ^  (4-  *i)  +■ 


(330 


and 


£2  —  ^  ^ 


(331 


Making  the  kappa  substitution  into  the  second  bracketed  term 
it  becomes 


[p,y  +m,  +h  +  zi~^j  W*  + A  -  kVj 


(332 


For  the  antisymmetric  cases,  •  '  Thus,  the  parenthe 

sised  term  becomes 


(£  ‘M  *  W  -2M)Xt  ~  2M'y)  Xz)  ( 333 


or 


(334 


Thus, 


+  +  M  +  TAKt-'O  Vi] 


(335 


<j=  ?>(*-»■)  +ih-e' 


(336) 


and 


£y—  A +  *rA  (337) 

Expressing  the  integral  in  terms  of  these  constraints 


yields 


X 

«w, 


f  £ 


(338) 


Collecting  terms 


T,  =  -&< 


—  +■  4*) 

W/  J 


X v 

L  ^  4-4,  f  ^  339) 


and 


,A 


I3=-v4«t  j ^  ^  ^  +(*>**  *■$*>/- &i)  ^ 


*-(<><>  +<; ty-c,^))?-  W*-*)*]  *fx  (340) 


m 


Integrating 


!,-*<¥(-  gg+ito  *<«-<*)£ 


+T fe<.i-4<V-44)X- 


(341) 


The  four  constants  can  be  expressed  in  terms  of  the 
nodal  values  of  the  fluxes.  The  first  constant  is 


A  =  —<5((-X'|  •+■  -J-  of2  [m3  —2m,) 


or,  utilizing  the  equivalence  between  and  j(3 


(342) 


dz  d 3 

A  (°(i  *2  +  "+  ®b) 

f  ~2  (*»*J  **  4j) 


(343) 


Thus,  from  the  finite  element  approximation,  it  becomes 


d,=  —7?+i-(a  +  4) 


(344) 


The  second  constant  becomes 


d2  =  «U  = 


(345) 


The  third  constant  is 


£3  —  W*  +  *  -** ~2Mj') 


(346) 


or,  utilizing  the  equivalence  between  y*  and  y3 

t3  =  -  -  b 


+ ~2  ( Pty  fb'Ut'-t  h) 

+\[h  +  h) 


Thus,  from  the  finite  element  approximation,  it 

<0  =  -ff  '+  T  M  V 

The  fourth  constant  becomes 

4=  +  i ft (V+^' )+ h  -  i- (</>,' *  4) 

For  Case  2 


and 


X?  =  ~  & 

Thus,  the  integral  becomes 


(347) 

becomes 

(348) 

(349) 

(350) 

(351) 


W  fat,  +t,t y  -  £,0*5  +  d,  C¥ 


(353) 


J3  =  -X*  ~j^T  l  )zc‘ 6  “&(**’*•  4 4  -  e*c7) 

-lofat)  -+  44  -  4 4>  4*4^ 


Collecting  terms,  yields 

J3=  ~~  W,  -lo&ctj  (354) 

For  Case  3 


and 


Xf>  **  X  -  *,  =  O 


—  xx  x,  ci 


Thus,  the  integral  becomes 


¥/*3' 

(,o*n 


yV 

+  is(c,ei+t,t9-tef)~~ 


+&&<>+ 44 - e,*)**- &ek6JX1}  )* 

^  o 


(355) 


(356) 


(357) 


or 

Ij  =  £—U 44+  c>  +  44  ~ 4 4) 

4  4+44  *44J  -  Jo 4  4  J 


(358) 
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Collecting  terms,  yields 


which  is  the  same  as  for  Case  2. 

Substituting  for  the  constants 

=  If  {?{-*<(!+<  +  4)(-i/;*4\ 4) 
-sHt+A.KtKti'-K# ) 

+s  (*4)(~2<P,'+  A'+  <?,') 
-io(i*)(rt  ■+<?;) } 


Rearranging  terms 

-5-  -r(4 *  4)(rt’+  4') 

ts(2t)(£\ <%)  i.  r(i/,)(4 +0 


Collecting  terms 


(362) 


-2  (fcrtXJ’W/)} 


in  matrix  form 


H  3 

3 

< 

3  / 

/ 

< 

.3  1 

1 

mi 

KJ 

(363) 


Therefore,  the  local  matrix  for  the  third  scatter  term,  for 
the  antisymmetric  cases,  is 

- 

3  3 
1  I  * 
i  I 


*  **' 
mL.  ==  <X  — ; 


3 
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